{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the main goals of BIDMat/BIDMach is to make model creation, customization and experimentation much easier. \n",
    "\n",
    "To that end is has reusable classes that cover the elements of Learning:\n",
    "\n",
    "* Model: The core class for a learning algorithm, and often the only one you need to implement.\n",
    "* DataSource: A source of data, like an in-memory matrix, a set of files (possibly on HDFS) or a data iterator (for Spark).\n",
    "* DataSink: A target for data such as predictions, like an in-memory matrix, a set of files, or an iterator. \n",
    "* Updaters: Update a model using minibatch update from a Model class. Includes SGD, ADAGRAD, Monte-Carlo updates, and multiplicative updates. \n",
    "* Mixins: Secondary Loss functions that are added to the global gradient. Includes L1 and L2 regularizers, cluster quality metrics, factor model metrics. \n",
    "* Learner: Combines the classes above and provides high-level control over the learning process: iterations, stop/start/resume\n",
    "\n",
    "When creating a new model, its often only necessary to creat a new model class. We recently needed a scalable SVD (Singular Value Decomposition) for some student projects. Lets walk through creating this from scratch. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalable SVD\n",
    "\n",
    "This model works like the previous example of in-memory SVD for a matrix M. The singular values of M are the eigenvalues of M M^T so we do subspace iteration: \n",
    "\n",
    "$$P = M M^T Q$$\n",
    "$$(Q,R) = QR(P)$$\n",
    "\n",
    "But now we want to deal with an M which is too big to fit in memory. In the minibatch context, we can write M as a horizontal concatenation of mini-batches (this assumes data samples are columns of M and features are rows):\n",
    "\n",
    "$$M = M_1 M_2 \\cdots M_n$$\n",
    "\n",
    "and then $$P = \\sum_{i=1}^n M_i M_i^T Q$$\n",
    "\n",
    "so we can compute $P$ by operating only on the minibatches $M_i$. We need to be able to fit $P$ and $Q$ in memory, their size is only $k~ |F|$ where $k$ is the SVD dimension and $F$ is the feature set. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining a new model class which extends BIDMach's Model class. It will always take an Options instance as an argument:\n",
    "\n",
    "<b>\n",
    "<code style=\"color:blue\">\n",
    "class SVD(opts:SVD.Opts = new SVD.Options) extends Model(opts)\n",
    "</code>\n",
    "</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The options are defined in the \"Object\" associated with the class. In Scala \"Object\" defines a singleton which holds all of the static methods of the class. It looks like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><code style=\"color:blue\">\n",
    "object SVD  {\n",
    "  trait Opts extends Model.Opts {\n",
    "    var deliciousness = 3\n",
    "  }\n",
    "  \n",
    "  class Options extends Opts {}\n",
    "  ...\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Truthfully, an SVD model doesnt need a \"deliciousness\" option, in fact it doesnt need any Options at all - or rather what it needs is inherited from its parent. But its there to show how options are created. The Opts are defined as a trait rather than a class so they can be mixed in with the Options of other learning classes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Local Variables and Initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three variables we need to keep track of:\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  var Q:Mat = null;                                        // (Left) Singular vectors\n",
    "  var SV:Mat = null;                                       // Singular values\n",
    "  var P:Mat = null;                                        // P (accumulator)\n",
    "</code></b>\n",
    "\n",
    "and an initialization routine sets these to appropriate values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minibatch Update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each update should update the stable model: Here its $P$:\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  def dobatch(mats:Array[Mat], ipass:Int, pos:Long):Unit = {\n",
    "    val M = mats(0);\n",
    "    P ~ P + (Q.t &ast; M &ast;&#94; M).t                               // Compute P = M &ast; M&#94;t &ast; Q efficiently\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Score Batches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The score method should return a floating point vector of scores for this minibatch.\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  def evalbatch(mat:Array[Mat], ipass:Int, pos:Long):FMat = {\n",
    "    SV ~ P ∙ Q;                                            // Estimate the singular values\n",
    "    val diff = (P / SV) - Q;                               // residual\n",
    "    row(-(math.sqrt(norm(diff) / diff.length)));           // return the norm of the residual\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update the Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of a pass over the data, we update $Q$. Not all models need this step, and minibatch algorithms typically dont have it. \n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    "  override def updatePass(ipass:Int) = {   \n",
    "    QRdecompt(P, Q, null);                                 // Basic subspace iteration\n",
    "    P.clear;                                               // Clear P for the next pass\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convenience Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're done defining the SVD model. We can run it now, but to make that easier we'll define a couple of convenience functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An in-memory Learner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><code style=\"color:blue\">\n",
    "  class MatOptions extends Learner.Options with SVD.Opts with MatSource.Opts with Batch.Opts\n",
    "  \n",
    "  def learner(mat:Mat):(Learner, MatOptions) = { \n",
    "    val opts = new MatOptions;\n",
    "    opts.batchSize = math.min(100000, mat.ncols/30 + 1)\n",
    "  \tval nn = new Learner(\n",
    "  \t    new MatSource(Array(mat), opts), \n",
    "  \t\t\tnew SVD(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\tnew Batch(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\topts)\n",
    "    (nn, opts)\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A File-based Learner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b><code style=\"color:blue\">\n",
    "  class FileOptions extends Learner.Options with SVD.Opts with FileSource.Opts with Batch.Opts\n",
    "  \n",
    "  def learner(fnames:String):(Learner, FileOptions) = { \n",
    "    val opts = new FileOptions;\n",
    "    opts.batchSize = 10000;\n",
    "    opts.fnames = List(FileSource.simpleEnum(fnames, 1, 0));\n",
    "    implicit val threads = threadPool(4);\n",
    "  \tval nn = new Learner(\n",
    "  \t    new FileSource(opts), \n",
    "  \t\t\tnew SVD(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\tnew Batch(opts), \n",
    "  \t\t\tnull,\n",
    "  \t\t\topts)\n",
    "    (nn, opts)\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Predictor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A predictor is a Learner which runs an existing model over a DataSource and outputs to a DataSink. For SVD, the predictor outputs the right singular vectors, which may be too large to fit in memory. Here's a memory-to-memory predictor:\n",
    "\n",
    "<b><code style=\"color:blue\">\n",
    " class PredOptions extends Learner.Options with SVD.Opts with MatSource.Opts with MatSink.Opts;\n",
    "  \n",
    "  // This function constructs a predictor from an existing model \n",
    "  def predictor(model:Model, mat1:Mat):(Learner, PredOptions) = {\n",
    "    val nopts = new PredOptions;\n",
    "    nopts.batchSize = math.min(10000, mat1.ncols/30 + 1)\n",
    "    nopts.dim = model.opts.dim;\n",
    "    val newmod = new SVD(nopts);\n",
    "    newmod.refresh = false\n",
    "    model.copyTo(newmod)\n",
    "    val nn = new Learner(\n",
    "        new MatSource(Array(mat1), nopts), \n",
    "        newmod, \n",
    "        null,\n",
    "        null,\n",
    "        new MatSink(nopts),\n",
    "        nopts)\n",
    "    (nn, nopts)\n",
    "  }\n",
    "</code></b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets try it out! First we initialize BIDMach as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 CUDA device found, CUDA version 7.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.9791425,12615962624,12884705280)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import BIDMat.{CMat,CSMat,DMat,Dict,IDict,FMat,FND,GDMat,GMat,GIMat,GLMat,GSDMat,GSMat,\n",
    "               HMat,IMat,Image,LMat,Mat,SMat,SBMat,SDMat}\n",
    "import BIDMat.MatFunctions._\n",
    "import BIDMat.SciFunctions._\n",
    "import BIDMat.Solvers._\n",
    "import BIDMat.JPlotting._\n",
    "import BIDMach.Learner\n",
    "import BIDMach.models.{FM,GLM,KMeans,KMeansw,ICA,LDA,LDAgibbs,Model,NMF,RandomForest,SFA,SVD}\n",
    "import BIDMach.datasources.{DataSource,MatSource,FileSource,SFileSource}\n",
    "import BIDMach.mixins.{CosineSim,Perplexity,Top,L1Regularizer,L2Regularizer}\n",
    "import BIDMach.updaters.{ADAGrad,Batch,BatchNorm,IncMult,IncNorm,Telescoping}\n",
    "import BIDMach.causal.{IPTW}\n",
    "\n",
    "Mat.checkMKL\n",
    "Mat.checkCUDA\n",
    "Mat.setInline\n",
    "if (Mat.hasCUDA > 0) GPUmem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll run on the MNIST 8M (8 millon images) digit data, which is a large dataset distributed over multiple files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "val dir=\"../data/MNIST8M/parts/\"\n",
    "val (nn, opts) = SVD.learner(dir+\"data%02d.fmat.lz4\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set some options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "opts.nend = 10;\n",
    "opts.dim = 20;\n",
    "opts.npasses = 2;\n",
    "opts.batchSize = 20000;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and release the beast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pass= 0\n",
      " 4.00%, ll=-0.01936, gf=4.010, secs=0.6, GB=0.13, MB/s=200.38, GPUmem=0.983898\n",
      "25.00%, ll=-0.00600, gf=13.319, secs=1.2, GB=0.82, MB/s=665.60, GPUmem=0.983898\n",
      "48.00%, ll=-0.00578, gf=16.846, secs=1.8, GB=1.51, MB/s=841.88, GPUmem=0.983898\n",
      "70.00%, ll=-0.00638, gf=18.716, secs=2.3, GB=2.20, MB/s=935.32, GPUmem=0.983898\n",
      "91.00%, ll=-0.00555, gf=19.860, secs=2.9, GB=2.89, MB/s=992.47, GPUmem=0.983898\n",
      "100.00%, ll=-0.00647, gf=20.827, secs=3.0, GB=3.14, MB/s=1040.82, GPUmem=0.983898\n",
      "pass= 1\n",
      " 4.00%, ll=-0.00599, gf=18.593, secs=3.5, GB=3.26, MB/s=929.19, GPUmem=0.983898\n",
      "25.00%, ll=-0.00425, gf=21.004, secs=3.8, GB=3.95, MB/s=1049.78, GPUmem=0.983898\n",
      "48.00%, ll=-0.00573, gf=23.088, secs=4.0, GB=4.64, MB/s=1153.97, GPUmem=0.983898\n",
      "70.00%, ll=-0.00617, gf=24.961, secs=4.3, GB=5.33, MB/s=1247.65, GPUmem=0.983898\n",
      "91.00%, ll=-0.00391, gf=26.602, secs=4.5, GB=6.02, MB/s=1329.75, GPUmem=0.983898\n",
      "100.00%, ll=-0.00589, gf=27.082, secs=4.6, GB=6.27, MB/s=1353.77, GPUmem=0.983898\n",
      "Time=4.6360 secs, gflops=27.07\n"
     ]
    }
   ],
   "source": [
    "nn.train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model matrices for this model hold the results. They are generic matrices, so we cast them to FMats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "val svals = FMat(nn.modelmats(1));\n",
    "val svecs = FMat(nn.modelmats(0));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGQCAYAAABYs5LGAAApA0lEQVR42u3de4xV9bnGcRNiGmOMIWka0zQnTZPGNE3/OGnCHydNYzxtaNqYXkzrDa0FUapyCipivFAq1aLHFg+IQksdVEqnWm5VsVZLq1gGLXW8U61QrB4VpXhBqMjlPbx7nRlnxhkYhpnZa+39+SY7m7X2dnh91o951rPe3/qtwwIAAFSew0gAAABDBwAADB0AADB0AADA0AEAYOgAin9Ehx3W71dv328mvv71rzft/zvA0AGGXnleffXVmDVrVlOfzAAMHWgwQ6cPQwcYOlABwwJDBxg60ASG3tf3d+3aFdOnT4+Pfexjcfjhh8cRRxwRxx13XO3Pff38vn7W/mrq+dmePXvi6quvjo9//OMxYsSIzu89+eSTcdZZZ8W//du/ddbz6U9/Oi6//PLYtm0bQwcYOsDQe/v+6NGjD/qy/WAY+uc+97kPfHfhwoU1c++rhjT2rVu3upoBMHSgcQ19IJPiFi1a9AGT3bJlS7z22mtDbug9X+3t7d22p06dWkvxc+bM+cB+hg4wdIChd9nXM53fd999AzLngfw3ixcv7nYJfcyYMd0+37lzZ21/mnrX/Z/4xCcYOsDQAYbedV/2p3v2tYfL0HsycuTIfv2/ZF+doQMMHWhYQx/I9wdqwENh6PvrnR+qITN0gKEDDW3oPRN6PQ39qKOO6vZ5zr6vl1YAGDpQKUMfNWpUt33PPvtsv35+z1vaHnrooVi3bt0hGfqJJ57Y7fNbb731A9/Zvn177ZY6hg4wdIChd/l+3n/edV/e/5199HvuuWe/P/9Tn/rUoN3q1sEjjzzS7fOPfOQj0draWkvqWdOqVatqt63115CtogcwdKBpDD3v6c4FXA5kel0XfkmuvPLKD3wnZ58fiqEneRvdhz70oUExYYYOMHSgaQw9eemll+KMM86II488snYp/fOf/3w8+OCD3b6bvfaeTJ48ubY/DXjs2LG1k4NDNfRk06ZNMWXKlFo7IH9+nkzk3/GZz3wmzj777FqSZ+gAQwfQD55++ulupjeQvjUAhg5gGMk13OfOnVtbHS5Zv359LaV3NfS8FA4ADB0o8z/GA1ySPvPMM4kEgKEDZefUU0+Nz372s7X7wLNXna+PfvSjcfLJJ8edd95JIAAMHQAAhg4AABg6AABg6AAAgKEDAMDQAQAAQwcAAAwdAAAwdAAAGDoAAGDoGHra29u7bY8fPz6+//3vl+r1ve99T020opW6aDXA16RJkxh6M/DHP/6x2/a3v/3t0tX4hz/8QU20opW6aDVA0tQZehMaep5dlo2///3vaqIVrdRFK4aOgzH0gz3wAIByw9AbnOydp5m3trZ22z9u3LjaJaOOs8x8r/f2smXLSlVPbre1tZWqno7trMvx6992x8vxq+Z4L+PxK+t4Z+hNmtD10KtbE61oRSt1SegMvRM99OrWRCta0UpdDJ2hD/jAAwDKDUOX0J3xSge0UhOtJHRU1dD10KtbE61oRSt1MXSGLqFLB7RSE60kdIbeSIauhw4AjQVDl9Cd8UoHtFITrSR0VNXQ9dCrWxOtaEUrdTF0hi6hSwe0UhOtJHSG3kiGrocOAI0FQ5fQnfFKB7RSE60kdFTV0PXQq1sTrWhFK3UxdIYuoUsHtFITrSR0ht5Ihq6HDgCNBUOX0J3xSge0UhOtJHRU1dD10KtbE61oRSt1MXSGLqFLB7RSE60kdIZeRdrb22tm3tra+gFDzzPMjkGZ77Zt27Ztu7rbDF1Cd8YrHdBKTbSS0FFVQ9dDr25NtKIVrdTF0Bm6hC4d0EpNtJLQGXojGbr70AGgsWDoErozXumAVmqilYSOqhq6Hnp1a6IVrWilLobO0CV06YBWaqKVhM7QG8nQ9dABoLFg6BK6M17pgFZqopWEjqoauh56dWuiFa1opS6GztAldOmAVmqilYTO0BvJ0PXQAaCxYOgSujNe6YBWaqKVhI6qGroeenVrohWtaKUuhs7QJXTpgFZqopWEztAbydD10AGgsWDoErozXumAVmqilYSOqhq6Hnp1a6IVrWilLobO0CV06YBWaqKVhM7QG8nQ9dABoLFg6BK6M17pgFZqopWEjqoauh56dWuiFa1opS6GztAldOmAVmqilYTO0BvJ0PXQAaCxYOgSujNe6YBWaqKVhI6qGroeenVrohWtaKUuhs7QJXTpgFZqopWEztAbydD10AGgsaisoW/YsCGOOeaYzu0VK1bE0UcfHYcffngce+yxte1m3C+hSwdqohWtJPTKGPrMmTNj1KhRcdhh75dz/PHHx7Jly2p/XrRoUZxwwglNub+/hq6HXt2aaEUrWqmrYQx9wYIFsWfPnm6GfuSRR9b2Jfl+1FFHNeV+CV06UBOtaCWhV8bQO4vpYugjRozo9lleim7G/T3JS/ELFy6M+fPnd9s/bty4aGlpiba2ttp2vtu2bdu27epuN4yhH3HEEb0m1mbbL6FLB2qiFa0k9Eob+ujRozsnht155521HnMz7u+voX9nzJjSDUb9O1rRSl20YuixZMmSWl85940cObJmcs24v7+GfslXvuKMVzqglZpoJaGjanzgPvT//E+iAEADwdAbnPb29pqZt7a2dtt/+b//e+2SUcdZZr7XeztvwytTPbndMemkLPV0bGddjl//tjtejl81x3sZj19ZxztDb9KEfsVnPhPxzjulqlH/jla0UhetJHQcpKFf/h//EfHMM6WqUf+OVrRSF60YOg7S0L+fK8vdcw9hAKBBYOhNaugXn3hixLx5znilA1qpiVYSOqpAX5PivrMvob91/vmlmtSRKx2VbZJJx3r5ZZuUk3WVbZJQGY9fb++OX7XGexmPX1nHO0Nv0oQ+aZ+Zx+mnR+zd64xXOqCVmmgloaOqhl478OeeG/HKK8QBgAaAoTepodfWcr/mmoi1a53xSge0UhOtJHRU1dBrz0P/5S8jbr+9NDW6B5ZWtFIXrRg6BpLQ16yJuO46Z7zSAa3URCsJHWWnr1nuaegPL10a702Y0Dk4yzjr1rZt27Ztm+WOAyX0fJ56Pkb13Xed8UoHtKIVrSR0VNHQaz30ZOrUiOeeK0WN+ne0opW6aMXQMZCEntx4Y8T99zvjlQ5oRStaSeiooqF3Hvi77or4+c8JBAAVh6E3e0J/8smIadOc8UoHtKIVrSR0VNHQO3vob70VceaZpahR/45WtFIXrRg6+qCv29bGjRv3/m0P48fHi489VvfbLvLBEGW7DaStra1U9XRsZ11lu22mjMcv3ztejl81x3sZj19ZxztDb9KE3u3A//CHEY8+SiQAqDAMvUkNvbOHntx6az7Pse416t/RilbqohVDx0EaemcPPXnggYjrr697jfp3tKKVumjF0HEoCX3TpogLLnDGKx3QilZqktBRNUPvduB37Yo47bTiHQBQSRi6hF6QCb3OZ5zSAa1opS5aMXQcpKF366En2UPv8Z3hRv+OVrRSF60YOg41oecs95zt7oxXOqAVrdQkoaN87O956F0XJnj1nnvinxdc4PnCtm3btu156Kh0Qt+ypbZinDNe6YBWtFKThI4KGfoHeuhJrumea7vXCf07WtFKXbRi6DjUhF6MhognnnDGKx3QilZqktBRFUPv9cDffHPxfHQAQOVg6BL6+9x/f8Tcuc54pQNa0UpNEjqqYui99tCfey7i4ovrVqP+Ha1opS5aMXQMRkJ/992IMWMi9uxxxisd0IpWapLQUQVD7/PAT5wY8dJLBAOAisHQJfTuXHddxJo1znilA1rRSk0SOspEXyvFnXTSSb2uNPTGT38a8ctf1mWlo5aWltKtvLQsl8SN8q0MlXWVbaWqMh6/3t4dv2qN9zIev7KOd4YuoXfn4YcjZs50xisd0IpWapLQUQVD7/PAv/pqxLnnEgwAKgZDl9C7s3dvxOmnR2zf7oxXOqAVrdQkoaPsht7rfegdXHppxPr1w16je2BpRSt10YqhY7ASejJvXsS99zrjlQ5oRSs1Segou6Hv98Dfc09EznYHAFQGhi6hf5Bnnom4/HJnvNIBrWilJgkdZTf0/fbQ33kn4owziglyw4j+Ha1opS5aMXQMZkJPJkyI2LxZOpAOaEUrNUnoKLOhH/DA/+hHEY88QjgAqAgMXULvncWLI379a+lAOqAVrdQkoaPMhr7fHnry0EMRP/7xsNaof0crWqmLVgwdfdDXw1nGjRu338X+X2pri+3jxw/rwwXywRBle9hB2z4dylRPx3bWVbaHVZTx+OV7x8vxq+Z4L+PxK+t4Z+hNmtAPeOB3744YMyZi507iAUAFYOhNaugH7KEnU6ZEPP/8sNWof0crWqmLVgwdB2noB+yhJzfcELFq1bDVqH9HK1qpi1YMHUOR0H/zm4iWFulAOqAVrdQkoaOsht6vA//44xHTpxMPACoAQ5fQ++bNNyPGjpUOpANa0UpNEjrKauj96qEnZ50VsXXrsNSof0crWqmLVgwdQ5HQkyuvjHjsMelAXbSilZokdJTR0Pt94BcujFixgoAAUHIYuoS+f/Ky0uzZ0oG6aEUrNUnoKKOh97uHvnFjxEUXDUuN+ne0opW6aMXQMVQJ/b33Ik47rVgKVjpQF61opSYJHeUy9IM68JMmRbzwAhEBoMQwdAn9wMyaFfHgg9KBumhFKzVJ6Cibofe7h54sXRqxaNGQ16h/RytaqYtWDB1DmdDXrYu4+mrpQF20opWaJHSUzdAP6sC/9lrEOecQEQBKDEOX0PtHXqLftk06UBetaEUrCR31oL29vWbmra2t3fafdNJJtR5Qx6DM9/1t/2vKlGhftKjf3x/IdktLy5D+/IFsL1u2rFT1dGxnXWWqp6zHr7d3x69a472Mx6+s452hS+j9Y8GCiJUrpQN10YpWtJLQUSZDP9gDH7/7XcRNNxESAEoKQ5fQ+8ezz0Zccol0oC5a0YpWEjrKZOgHdR96smNHxJgxEXv3DlmN7oGlFa3URSuGjqFO6Ml550W8/LJ0oC5a0YpWEjrKYugH3UNPrr02Yu1aYgJACWHoEnr/yVvffvUr6UBdtKIVrSR0lMXQD7qHnrS1FSl9iNC/oxWt1EUrho7hSOjZPz//fOlAXbSiFa0kdJTF0AfUQ88Z7jnT/V//IigAMHRUNqEneS963pMuHaiLVrSiFUNH/Q19QD30JFeLy1XjhgD9O1rRSl20qoihjxo1Kt566y3uWuWEfvfdxbru0oG6aEUrWjWvoT/++OPx4Q9/eOBmgkEz9AH10JOnnoqYNo2gANDMht7BwoUL4/DDD49FixY5AlVL6G+/HXHmmdKBumhFK1ox9Pf55Cc/GSNGjOh8pcljeAx9wD305OyzI15/fdBr1L+jFa3URauKGfpNN90UH/rQh2LZsmWctmoJPbnqqoh166QDddGKVrRqVkNfu3ZtjBw5MqZOncph62zoA+6hJ9kqWbqUqADQrIZ+3HHHxbZt26he9YT+4IMRs2ZJB+qiFa1o1ayGjvIY+iH10F94IWLSpEGvUf+OVrRSF60YOoYzoe/eHXHaaRHvvScdqItWtKIVQ0c9Df2QeujJhRdGbNxIWABg6KhsQk9mz84fKh2oi1a0ohVDRz0N/ZB66Mny5RG33DKoNerf0YpW6qIVQ8dwJ/T29ogrr5QO1EUrWtGKoaOehn7IPfStWyPOOouwAMDQUemEnowdG/Hmm9KBumhFK1oxdNTL0A+5h55Mn56P0Bu0GvXvaEUrddGKoaMeCf3mmyN+8xvpQF20ohWtGPrgs2HDhjjmmGM6t1esWBFHH3107Uluxx57bG27kff319APuYee/P73ETfc4GwJABj64DJz5swYNWpUHHbY+2Uef/zxnU91y+evn3DCCQ29f1gT+vPPR1x8sXSgLlrRilYMfXBZsGBB7Nmzp5uhH3nkkbV9Sb4fddRRDb2/J5ncFy5cGPPnz++2/+tf/3q0tLREW1tbbTvfD3b7ttT71FNrS8EO5L/vuT1r1qxD+u+HYvu2224rVT0d21lXmeop6/HL9+x1On7VHe9lPH5lHe8N2UPvaugjRozo9lleom7k/cOa0JPx4yNWr5YO1EUrWtFKQh9aQz/iiCN6TbKNur+/hj4oPfTkxhsjzj034pJLIh54IGLXrgAAMPRBN/TRo0d3Thi78847a73nRt4/7Ak92bs34i9/ifjhD4vE/qtfRbzxhnQgSamLVrRi6INn6EuWLKn1m3PfyJEja+bXyPv7a+iDch96b7z0Uk5miDjzzOIBLn/7W7//U/fABq1opS5aMXT0Tnt7e83MW1tbu+0fN25cbUB2nGXm+2Buv7B+ffzt+utj14QJEZdeGq8tXRp//P3v9/vf56z9oapnoNsdk07KUk/HdsdkoTLpVcbjl+8dL8evmuO9jMevrOOdoTcJQ9ZDPxB5OX7duuJBLmefHXHHHYO6XCwAQEJvakMf1B56f3nxxYif/rS4HD9nTnEfexf076pdF61oRSuGjjoY+pD10PvDO+8US8bm7PjLL4/4059q97Lr31W7LlrRilYMHc2S0HuSt9098kjxkJdzzontORiffbbY7yxckqKVumjF0HFgQx+2Hnp/2bgxYurUiAsuyMsHEVdfncvcFZflS2TwAFBWGHqDU69Z7oc0a/Ttt2Pzb34TL82YETsnTqwZ/I4rrogNP/lJvPzgg7XL82a5m+VulrtZ7sZ7922G3qQJva499D7osye1z+Dj4YeLx7VedFHEGWdEXHVVxPLlEc89VzP4Ya+prFqpiVa0atq6GHqTGnopeug96HdPatu2ovfe0vK+wecqdfn0uUE2eP07WtFKTVWpi6E3qaGXrod+KOSs+TT4hQsjpkyJGDOmmEG/Zk3Ee+85+ACaAoYuoTfeGe/WrUVaz8Vsxo4tlqLtcc+7dCBJqYlWEjoawtAr1UM/FF5/PeLXv444//xiBn3e/34QK9Xp3wWtaKWmitTF0BucSs5yH6Kf/+gvfhHbrr22tlLd9mnT4qmf/zw2/X9yN8vdLHez3M1yN8sdlUzoDdVDP1jefbd4dnsuaJOX5HP2fN4HDwAVhqE3qaE3dA/9YHjtteKBMXlJ/sIL8+Hy3S7J699VuyZa0aqZ6mLoTWroTdND7y/5VLinn46YO7d4eMw119Tufc9HvpYRvU5a0UpdDJ2hS+gHIi/Jp17Tp8eevAVuxoyIJ54Y0gVsJClaqYtWDB0DMvSm7qEfDOvXRyxaFHHppUW//YYbinved+6kDQCGDgm9kungn/+M+O1vi/vbs2Vx3XURubb89u20kqRoRSuGjqGlr9vWTjrppNLdBtLS0lK620Dy1pTePn/h6adj/bx5xSNfzzgjdlx2WTw7d27848knh6W+rMvx6992z/cy1FfG47e/8e74VWO8M3QJ3RnvodaUPfe8f/f664sJdVdcEXHXXcUMelqpi1a0ktAxlIauhz5E7NoV8eijEfvSe5x1VsTFFxcr1f3jH7QBwNAhoVeyprwV7plniofG5MNisu++eHFh+rQyrmhFK4aOwTB096HXoabVqyOuvjriu9+NuO++Q7oNzv3CtKKVuhg6Q5fQ613T3/5WPL/9vPMiVq2K2LOHVsaVumjF0DEwQ9dDLwF//Wuxnvx//Vdx+1teogcAhg4JvaI1PfVUMTN+0qSINWv6ZeySlHFFK3Ux9CbDfeiHtj2c9+U+fuut8W4+IOaii2LzXXe5D9196O5Ddx/6QW0zdAndGW/ZavrLX4rb3fKVf6aVcUUrdUno6MvQ9dArwMMPF490zXXkH3+cHgAYOiT0ytaU/fTsq2d/Pfvs2W+nlXFFK3UxdIbegfvQK1ZTGnvOhJ84MeIHP4hn5s8flAVqGv34GVe0aqa6GLqE7oy3SjXlYjT33Rd7Tz014rTTIi64IGLWrIglSyL+/OeIzZvrdvubJEUrWjF01MHQ9dAbgEzomzYVyT2f2Z6r0E2YUHv6W63vnuvJr1xZXKZ/+216AQ0OQ5fQnfE2WjrI57OvXx9x770RP/tZxLRpxVPgxo+PmDEj4pZb8nphxIYNETt3SlJqopWEjiobuh56dWsacF1bthRPglu+PGL27Nr97nHKKUWinzkz4tZbI+6/v3igzJtvNrdWxpWaKlgXQ5fQnfE2czrIhJ5p/pFHCqO/8caIyy+PGDu2SPV56f6GGyKWLo1Yu7Z4DGwfk/EkKeOKVgwddTB0PXQckOy753rz+QCZ7NFfe21x+1xOyDv//KJnn5fv88lxTz8d8frrNAMYOoaKvpZ+HTduXOmWUsxlJ8u2lGJbW1up6unYzrrq9vfv3h0v/fnP8URLS2zNZ73fdFP8a8qU2HXKKbEnL9/v2/fKvpOAP5ZEr46X41fN8V7G41fW31cMvUkTuh56dWsqtVYvvBDx618Xy9buO2lMs4916yLee49WxhWtJHQMhaHroVe3pspo9dprEXffnb9l8gwy4ic/iXjooWIWPq2MK1oxdAyOoeuhY1jJfnz24nM2fV6W/+EPI373u4g33qANwNAhoUsHlazp3XeziVzcOpcz6XNW/YoVEa+8Qit10Yqh42ANXQ+9ujU1lFa5lO1jj0X89KcRZ59dLGWbEzifeIJW6qIVQ4eELh1UsqZcg/6554pb5E4/vbgXPlN8jt2tW2llXKmLoaM3Q9dDR+nJSXW5cl1Opktznzw5oqUl4i9/KS7bA2DoDF1Clw4qVlOm9+efL1asmz69mFiXa9TnLXKZ6vfsoZVxJaEz9OY0dD306tZEqygSevbec/35XJM+J9ddd10xc/7VV2llXDVlXQxdQnfGKx1Uv6Z8mEw+RjbXnc/JdeedFzF/fm02/Yv5QBpaGVcSOkNvVEPXQ0dD8+KLxaI2V10VcdJJET/4QdGP91x4NDAMXUJ3xisdNLRWm559tnhS3KxZxYp1uajN738fsW0brYwrCR3Voa+Hs5y0L7WU7WEHLS0tpXvYQT6AoUz1dGxnXY5f/7a7vtfMva0t3pkxI3addlrsyEVtVq2KF555pumPX1nHe893493DWSR0CV06oFV3cmLdmjXFZLpM7vk42Px3MgxrzRtXtJLQMWiGrocOdOFf/yoeHPPf/12Ye645/8ADETt20AaVgaFL6M54pQNadSVNfPXqiGuvLcz9mmuKGfSDaO7GFa0YOgbN0N2HXt2aaDWMdaWJp5lnYs/laM8/P2L58oi//vWQnvFuXNGKoUNCVxOt6lXXli0Rd90VsXBhxCWXRIwZUzwpLteeX7fuoGbNG1fGFUPHoBm6HjpwiOSkuiefjLjjjogZM4rlaHO9+VzQJvvvmzfTCAwdErqaaFW5unI9+Y0bI1auLB4oc845xap1+efct2FD55rzxpVxxdAxaIauh17dmmhVobryiXHZg8/nveez3jPFX3llvJy/eHMluzwBKNHqdcZVteti6BK6M17pgFbDxTvv1B7/uj2fGHfZZRFTpkR85ztFP37SpGIVu7xkn0+VyxOBZ54pTgp27zau1MXQ0buh66EDJSLvg//HPyLyQTL5xLjFiyP+53+KSXcTJkScckrEd78bccUVEbNnF5/n99rbiz5+niig6WHoErozXumAVmWvKRN6JvVM7DnhbsmSiHnzikSfl/FPO614hGxe1s8V7zLl57Pi8999Gv7LL0fs3GlcSegMvRENXQ+9ujXRila9kil906baJf3OlD9nTkRe3p84MeLUUyPGji0u8+d99QsW5KLy71/a37w5Vud/5/hVti6GLqE745UOaNUsNb31VjER75FHIu65p7iHPi/tT5tWXNL/5jcjLr444mc/K9J9JnvHT0JHuQ1dDx3AB8jV7557rpiBf/31EeeeW0zay8v4eb/9Y4/p15cYhi6hO+OVDmilpr7revPNItH/4hfF5fvs2eeM/BtvjLjvvuIy///fX29cMXTUwdD10KtbE61oVde60rzTxLPfPnduYe5p8mn22bf/85+LkwBaMXRI6GqiFa0qVldehs/L8bffHnHVVe/PvM9Z92n0P/5xsbhOa2txOT8n4uUtd88/X8zez9v2jCuGjoEZuh46gCFj166Il16KeOGFiKeeilizJuLee4tb6Vpainvp85a7qVOLPn0urJOz8HOp3IsuivjBD4olc3Mmfs7U/9Wvisl8+XPB0JuV9n1nwWnmrXlm3IVx48bVLhl1nGXme723ly1bVqp6crutra1U9XRsZ12OX/+2O16OX7nH+wP33x//yJS/aVO88rvfxdM/+1ls+cUv4p0f/Sh2XHZZvHPOObE3Tf/CC2Pbvn3Pz5pV+14unWu8F9sMvUkTuh56dWuiFa2aVquchZ9JfdWqIumngeWCOpnscyZ+9vDzasD//m/E3r1NpxVDb1JD10Ovbk20ohWtepB9+JyMl7fWXXddxPnnR5x+esSllxa9+7zc/9e/Rmzf3tBaMfQmNXQ9dAANzY4dhYn/9rfFMrmXXBLxrW8Vq+Wl0c+aVSysk5/n6nq5lv5+JudVAYYuoUsHkhSt1NQcWuWa+P/8Z2H0Dz1ULH2bCT4v10+eXEzOy4V0crW8a68tLuvfdVfEww9/4FG3fdaVl/rzxOCNNyJeeSW/GLF+ffHgnWwHZLtg5criiXrZIrj55uKkgqFjoIauh17dmmhFK1oNIWnaGzZErF1bmHkabpp7roGf/fq8lL/P/HfmrXn59Lvcn2vljx9fnBDkVYD8Tm7n/vw8v5e38+Xte7kgT/7MNPM09TT3nAzI0CGhS1JqohWthpG8lL+vpq0//3mRvLO+TOKZyDOZD/FkPIYOPXQAaHAYuoTuTFySopWaaNUAdTH0JjV0PfTq1kQrWtFKXQydoUvo0gGt1EQrCZ2hN5Kh66EDQGPB0CV0Z7zSAa3URCsJHVU1dD306tZEK1rRSl0MnaFL6NIBrdREKwmdoTeSoeuhA0BjwdAldGe80gGt1EQrCR1VNXQ99OrWRCta0UpdDJ2hS+jSAa3URCsJnaE3kqHroQNAY8HQJXRnvNIBrdREKwkdVTV0PfTq1kQrWtFKXQydoUvo0gGt1EQrCZ2hN5Kh66EDQGPB0CV0Z7zSAa3URCsJHVU1dD306tZEK1rRSl0MnaFL6NIBrdREKwmdoTeSoeuhA0BjwdAldGe80gGt1EQrCb06bNiwIY455pjO7RUrVsTRRx8dhx9+eBx77LG17UbY319D10Ovbk20ohWt1NW0hj5z5swYNWpUHHbY++Uff/zxsWzZstqfFy1aFCeccEJD7JfQpQM10YpWEnrDGvqCBQtiz5493Qz9yCOPrO1L8v2oo45qiP39NXQ9dABoLJqqh97V0EeMGNHts7x03Qj7e9Le3l4z89bW1m77x40bV7tk1HGWme/13s4rDmWqJ7fb2tpKVU/Hdtbl+PVvu+Pl+FVzvJfx+JV1vDetoR9xxBG9Jtyq7+9vQtdDr25NtKIVrdQloXcx9NGjR3dOJLvzzjtrPelG2N9fQ9dDr25NtKIVrdTF0LsY+pIlS2p96Nw3cuTImik2wv7+GroeOgA0Fu5DbxIkdOmAVmqilYTO0BvQ0PXQq1sTrWhFK3UxdIYuoUsHtFITrSR0ht5Ihq6HDgCNBUOX0J3xSge0UhOtJHRU1dD10KtbE61oRSt1MXSGLqFLB7RSE60kdIbeSIauhw4AjQVDl9Cd8UoHtFITrSR0VNXQ9dCrWxOtaEUrdTF0hi6hSwe0UhOtJHSG3kiGrocOAI0FQ5fQnfFKB7RSE60kdFTV0PXQq1sTrWhFK3UxdIYuoUsHtFITrSR0ht5Ihq6HDgCNBUOX0J3xSge0UhOtJHRU1dD10KtbE61oRSt1MXSGLqFLB7RSE60kdIZeZR577LFu2xdddFHt4Ht5eXl5NcZr0qRJDB0AgGaDoQMAwNABAABDx6CxYsWK2kS5Mr1aW1vVRCtaqYtWA3z1nCvF0JuEnrPe1VSdmmhFK1qpS0JHJwd7JqcmddGKVmoqf10MHQCAJoOhAwDA0AEAAEMHAAAMvdnZu3dvzJ49O2bMmBHTp0+vvd5444261rR06dLacoVf/OIX40tf+lKtttdff700mt17771x3HHHlaaejRs3xtSpU2P06NG1V73JY3XFFVd01jNlypR4+eWX61bPpk2b4stf/nLpxn1vddV77PdWUxnGfl911XPs91ZT2cY+Q28yFi1aFLfffnvndt5Hef3119e1psmTJ8e6detqv3R3794dt912W0yYMKEUeuUs1ksvvbQ0hv7qq6/GmWeeGY8//nhpxtTYsWNraxrk8cvXkiVL6vYkvzxOHa8yjfu+6qrn2O+rpnqP/b7qqufY76umMo19ht6EnH322bFly5ZuZ5innHJK6er8whe+UPcaXnzxxVp62rFjR2kMPU3oT3/6U6mOVabL/GVWpuPX83iVZdz3ZxwNt3a91VSGsd/z7y3D2O9ZUxnHPkNvInpepsrBWIbLtl159NFHa8mlnuTl2HPPPTe2bt3a71/Ew0GaUP5iy8uzedwuu+yyUrRMbrnllli5cmVs3rw5Vq1aFddcc02pfvGWZdwfaBzVY+z3rKksY7/n31uGsd+zpjKOfYbeRPT2j7NM/eHVq1fH+PHj695Dv/DCC2v9urJpdPzxx3e7RJu/TC655JK61pS/yCZOnFj7RfbNb36z1lN8++23SzXOyzLu9/d31mvs96ypLGO/599bhrHfs6Yyjn2GLqGXorabb765FImz4x9ub68yGHrP45eX/epJprkNGzZ0bmev+mCfydzsCb2eY7+3k58yjP3eDL3eY79nTWUc+wy9ieh6KS3JNPCtb32r7nXlBJw8+67SlY16kBNuuvaCk6985St1ram3nmHZeuhlGfe9jaN6j/0Dje2yJPQyjP2eNZVx7DP0JmLx4sW1Gb4d5J9nzZpV15rWrl1bilReBUPP4zVnzpxOvbLneuONN9a1prxVJy89dqSmjrRZpuNVlnHfs64yjP2qGHoZxn7Pmso49hl6E1HG+9DLenm7jIaezJs3rzYxKC835i+4Xbt21bWeHD85jrKefGU/sV59xL7GUb3HfV911XPs9/fvredtaz3rqtfY76umMo19hg4AQBPD0AEAYOgAAIChAwAAhg4AABg6AAAMHQAAMHQAAMDQAQAAQwdQeu6///746le/WltFrDfy8Zv5eX4PAEMHUGJyedBcmvOOO+7otj/Xc8/98+fPJxLA0AFUgVxfO817zZo1te1Vq1bVtnM/AIYOoCLkQzsmTpxYe+b58uXLa+/nnXde7Ny5kzgAQwdQJd5666045ZRTasn85JNPLv0jegGGDgC9kAbO0AGGDqDCdFxyz2dU5+Q4l9wBhg6ggkybNq2WzP/whz/UtlevXl3bzv0AGDqACtBx21reptaVTOq5Pz8HwNABlJiOhWVyAZneyAVnLCwDMHQAABg6AABg6AAAgKEDAACGDgAAQwcAAAwdAAAwdAAAwNABAGDoAACAoQMAAIYOAAAYOlCFf4yHHVbZFwCGDqCLoasbAEMHGLq6AYYOgDEydIChA2DoABg6AIYOgKEDDH2AHHfccQwdYOgAmsXQ9/c9hg4wdAAMHQBDBxh6f8146dKlcfLJJ8fo0aPjsssui7fffrubUe/evTvmzJkTX/rSl2qv2bNn1/Z1fKfri6EDDB1Afw19+vSBv3ox9MWLF8f27dtj7969NXOfNWtWN0OfP39+3HHHHbXP87V8+fKYN2+ehA4wdACHZOhPPTXwVy+G3pVM3l/+8pe7fXbiiSd2JvJk165d8Y1vfIOhAwwdwCEZ+iDSmxl/4Qtf6PbZ/r7D0AGGDqCEhp5J/Gtf+1q3zzKN90zoHd9h6ABDB1ASQ9+4cWPtz9kfv/vuu2Pu3LndjDp76K2trZ099Oynd+2hn3DCCfHiiy8ydIChA6inoU+ePDm++MUv1maw52z2nTt3djP0TOc5sz2/k6/8c6b0DlauXNk5A56hAwwdQJ0MvYp1A2DoAEPvQiZuhg4wdAAVN3R1AwwdAGNk6ABDB8DQATB0oAkNvaovAAwdAAAMAv8H/+MRQ58tKJEAAAAASUVORK5CYII=",
      "text/plain": [
       "BufferedImage@70f39af4: type = 2 DirectColorModel: rmask=ff0000 gmask=ff00 bmask=ff amask=ff000000 IntegerInterleavedRaster: width = 500 height = 400 #Bands = 4 xOff = 0 yOff = 0 dataOffset[0] 0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "semilogy(svals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see how well we did, we will compute the SVD directly by computing $M M^T$ and computing its eigendecomposition. Normally we can't do this because $MM^T$ is nfeats x nfeats, but for this dataset nfeats is only 784. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..........\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "14.375"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tic\n",
    "val MMt = zeros(784,784);\n",
    "for (i <- 0 until opts.nend) {\n",
    "val Mi = loadFMat(dir+\"data%02d.fmat.lz4\" format i);\n",
    "MMt ~ Mi *^ Mi;\n",
    "print(\".\");\n",
    "} \n",
    "println;\n",
    "toc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we call an eigenvalue routine to compute the eigenvalues and eigenvectors of $MM^T$, which are the singular values and left singular vectors of $M$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "val (eval, evecs) = feig(MMt);\n",
    "val topvecs = evecs(?, 783 to 784-opts.dim by -1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenvectors have a sign ambiguity, and its common to see V or -V. So next we compute dot products between the two sets of vectors and flip signs if a dot product is negative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "val dots = svecs ∙ topvecs;\n",
    "svecs ~ svecs ∘ (2*(dots>0) - 1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets now look at the eigenvectors as small images, decreasing in strength from left to right. First the reference eigenvectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAACoCAAAAAD/IJEYAAA5TUlEQVR42u3d+bNuR13v8fvP6A9WaWkFsCRkHk/Gc3IIJAFCAgQCCBiGe5kTFbxeNWG8DGEOhCEQQuYQEsIcLSnUUkur/MF/5b7u867Td7meZ++zz0m4tbW6f3hq7bV6ffs7fXr3t9e3u//bHf/Jy+1b5aT1x1vPS9MHZ+zgTJ4GwQO2snq0fX0abG/TXFF7LpLupeptxvbh9pTkPT1n+G8TSBNIE0gTSBNIE0gTSM8rfg6ii30Msw2PfSyxbb/T4HD7rYO77GkD6SCy7++jO1k96aNT7XT2ornTTDul3gfe2xX2ssjB+9wJpAmkCaQJpAmkCaT/2kDax8mel7Hpc4lA9rLKwQmeXlhyegR3usv+0cLzEuOdhrn3Uu9px3inFL8dUIH7M3lKDC9rTiBNIE0gTSBNIE0g/RcG0v5+c9rhzV7e/8FNObg/nSpXz+PM+2lEU/uEZwf3iQMq4Tna+vTc44OLckrdzV5vPfcw8jQcYwJpAmkCaQJpAmkC6T81kP54UU5vFnWlwQ984AMrvewzCl+20ruV1eux9yebsmR1fyBtm2cl8pLac/8ivq2Knf50QA9YMrks22Y6KYWDPH2OOQHDfO973/ve+973+n3/+9+/coC9PMdb7z9RPrB32QmzvUTbdpWDB7oTSBNIE0gTSBNIE0iHHEh/vKvU0p/+x7Jtp/1HpUMLqe8973nPu9/97v+xKS7cGdo8iBkGncyQWt1XAWMf/vCH/2xTXGB15f07URRvSzOkDa9/6EMfiuCgdhA3+uM9ymi6FlOFsvKnvTIGVrBZsvc/N+XPTxR3PNoJp72caelS26YfpA4IpCXZKHxoU1x4SnxS//dN4QDEX1p/p/+krqXb7CwpE/2T8pl9P7wpQ1dDxtHXLx1sAmkCaQJpAmkCaQLpPwuQdgImyYdfZie/Kx89KZBwQLYU9453vOO2225761vf+ocnytve9ra3v/3tVKPOdgyzYjIbUKj6qKVQlN1xXwW8/eVf/uVHNuXOO+/EcJbbC0gN1gP20gZIkd3rqN25KS7Ijtpy8n0vTQ7trcqAooa0+M53vvO2TaGWNLDPxG7dxLAI3v7iL/7ir/7qr/B21113feRE+ehHP+rXHU9jeB8sbXMbZpYYqEW/KufQJ+314nP4THxWXLuPFJPlDFkfzZ1AGgHVUNcfnSipzuu3nSiu1Vkic6+M4VCEvf+1KSlqeNeAa262RPsE0gTSBNIE0gTSBNLhB1La3MZMbZCcD6WFXMqF+6otjbSX+eMMK+QEm1tvvfV1r3vdzTff/OpNuemmm/z5xje+kV4wjdHtmdYlk35pPGpv25QUSn1g4CluP/GJT9x9992f//znP/vZz37sYx8jSx6wsv1AePBe2sAj2kCKU37yk5/81Kc+5de1O0TOrn738SGNprfKX2xKBqNk/HidvJh/05ve9IY3vIFadC7vete79sESynSODhNg5uMf/ziu/vemuPjEieLPT22K67qSMLBX15nHryKr0VmMQCI/w16d19LvtynXB2md/vH26U9/mi0Y5TOf+Yw/gVxDNEDhb3nLW3gF5XOSbSCFovxHHZW5Sup685vf7M+3nihvOVGoNGfYBtLKTKPPxSe7qK8VHqWJ6CMbYnMz1iF176bPCaQJpAmkCaQJpAmkQwikoanGi2FmDLsZ7GOLwnjM49d9NbGC6M5x7Zg9pIX0RX74eeUrX3n99de/fFOuu+66G2644cYbb3zta1/Ln1RTv3H5mIVcMemXLSGHtLfccsvrNuX1r3+914FBfZwz2L333nvfffd985vf/MIXvuBOjlvwswQSvXALanrTplAcHjCsPj8gL/N/6Utf+vKXvwyWXBMpDGRXvyNsW6GIGqkun87FkcpgdEtvKqhGP2yjXWp51ate9ZrXvIaWihZ2JjrQNqOgWTehuOCjzOF+PV3t8tcq4Nkdlh1m2u46w3yMVZOGMVbPkkfWu+Vq9Vzu7wMkrGoXG7T3ta997eub8tWvfvWLX/yim7ShRdVQ46+8lgkQHPpc9cKY0aiaXChzu8bDuzaF3fP1EOXCndjeJ2+YHZmY3uof6RxNznns2LGrrrrqmmuuYRE+FsHQrkXMLP1zAmkCaQJpAmkCaQLpsAEphdImopmBkGOc/ekTpTF3Nz+zKS7iYIy/V3FRyQfkCUWcHmDABngGhPAtTOJGhUkhoS/NGX7I34iW/3mEJudD7bpNecUrXuEaBa1olGoY7zvf+c6jjz762GOPfetb3yIOCk1tL6cvm4OmnRDOjweYtYIOwUHoG5tyzz338EsuS96cyW9xy8BSsQEm6ZCLDwTm602gj7DBHQpnEi1q+mWb4iIfjc+lStVUH7C5Y075uc99jgm0SOo+6hMHP4zCUlxWHQxo3YuqbUcIY2K6ThN7dRNchz6Frzo47jX8kmZcsxRFJf7ORM8iVe1i4Nvf/vb999+vUwMnnRq3IT4l1JVoi6cWH2K+jy4jd6QIFlTwgw2N+nVNUtXq8SFfTXdUe+eJgrdUEbWdWdGeerfpeArkn/Dzkpe85IwzznjRi1504YUX8k9clXYTkgNSIWVh5ATSBNIE0gTSBNIE0iEE0pipZFpyNk3ZCPsLm9JY3H348dv9xt8MMMbfq/xCqjfcJxVuDGpZ5ZWbAjnM4w57jBlMvxSkPlJUPGZsR/Ioz27alwyk4vfHjx8/evTotddeiyb/QzBjYIn3PPDAA0899dQzzzzz0EMP4bOZVvSXc9YZAA9GwwCJMWygUFtFI0gJtPhBWCI+eXMmmi2iG87Ul3IeycX5Oh/yS3VINWmuRaIhS4FhDKI84iIvfelLr7zySmy4HnwOfbpG2YvcUR+BGbagIjyQnQaIQKvG8XwIS3jQRF2AC38ugTT0OcJgJd0yHEMwFvVeccUVV199tX6KbjXUVwH0+b2uhxXqhVcoKr8ENe3qxVIChWioef9crq8shUkMJx7zbhkzA0hIMVBfTfJpL2IVZTEwsroqCqSZUmDfe6KQIkdSqHd7uVd/FiW6pnNB0dlnn/17v/d7v/3bv/2CF7zgyJEjZEyQJusL9VGuA5pAmkCaQJpAmkCaQDqEQGpO2QMjRYEEL8FiI/tR/Gm0HZYUGvHnVzbFNezF6/aUZZ+im1lmCWZeZjA0wzumO1tq4g6lFEsMhY7hLEtojgw0a/zaHCXv77t4c7WaLkb67ne/+4Mf/OBHP/rRww8/TNeYARIA4A1D9twCV9duigvGIE7pnhhwDYQG+t86UVx7ig73JUu5FImQd9Kk2ICNv/e97+GBPjHcHLf6OCdLNOGTl1Mjz8A8IF166aXHjh0zWB9f0IfVEacTCseDF1HQNPBcfvnlRvOXXHKJF5uuzb/HNwCV2SskL8ODJvebnQ9FBRteR+T6669HU8Bw7rnngjegUjI2Suzo2z8jFtJs5xOrGZDAXusa8ojS6jeZaRi6zwlhoFSS+uVl0knBCQr6Dp3I448//sQTT+gf+7aBOLXXcQ+nQiqXDr1Fd0v/LHqvPkFgRtyu47jooovOO+88vxxMHS3SISlcoFM+Rxk/E0gTSBNIE0gTSBNIhxNIY7q2T/hf3RTWDUIDVGGJIf36E2kNeMoSKXc5ZVm6QCjqk3Bj3BITm8ZdJly6DjlBml39LoHklaZ0yeMCWWKHIlZp7p5gRSyFBxzukUce+f73v0/1iDf/7sWRjtDg3uuIIEVxrEUVuEKhRxSnOdTY7P7772c82ODEBOFqAgl1lhP0nIAmg7HKtOQOnxP5ACpHx4ZqpCiCCpnERweQLr74YjEJPDcWXw7oMUwn9C9Oowfupb4x/e///u+fddZZl112mSZG3FJfoBpjYSalLb8l0CcPdr/5aNd0WEIqCowlLrrgggvOOOOMF77whbiiHCKrFjV6LhFhrBsfHZOb5d1qkUdhOBQ1ny6yjUkvup87cfdmk7HddPYKSLSBOK/TLf7qV7/62c9+JgD27lhwlen7tpEV+nShDmsGs6HPyPJPbJTC4hVP+0iju+dadItPzLAR3OoTKQrPiVYa9wTSBNIE0gTSBNIE0iEEUp+3KYg2R3Ih07pu4tsFUDGeX1hyp+/lzau66UVuOqaqY3R8Caa+Mgv76N5A3yte9HrDfS329b2v3X77xL4EEvpUk1uogDJL80vWJW0IB79mVFMiPlG+7777XCBCQfyepsqMHN+8DdzROX78uJgH52Fb0zjvo75rUpcnIeh67LHHeJL7PP5lL3uZaiMRNrSPAEnTOIEisZx44+jRowGPjE1Mk73sBG/pZRgPHvzeeuutvGeM6euY/Oos+BwxicDY8POCF7zgzDPPPHLkCFRriJKJ3zolWmo1FB7KVyhnd8Ry6lB7sPRnmSjFEuyFVUD63d/9XUASOWBJ6yR6+umnOTHKpfYugVSA1DoxnBTRlQtC230wQIfJvKLRkh74gDrNvJep0Jz48iOKa5j8yU9+8m//9m9///d/7y11WnA+8gxaHt+cPv/hY/kqjZWr0dx9AoZ2rHrE4v2TwKfX8aljrVdlaEbXFxOcVmGPmZBqgdkE0gTSBNIE0gTSBNLhBJIBpUp333138OAElIuJZrqLmpoBLHDy26Sqyv7cjpHGAqSWfzcMbYG0YTGOR/Si8MsHH3wQKc01A9uCqOXWP34zPB7UpBpaMOA2lHefFsiJJborUmqZNyIlarAN4/XtfwmktsVBiuL4JRv70yv47GZRh2pExmST6TjnmjfeeCM/gyWIykJFm006Uw5ju2Aq0df555/P1zHA8NrFfzqkvboh7qL14igeXPjRp/fhTIV2Lc+G8EsvvRSEzjnnnCuvvDK5VCA+t2i1j4ZK4oy9EogHkNDRNFtjg8gqtMy7MFVl6gX+P/iDP7joootc06eu5Je//OU///M/0wMxW6fUSvUBJE1Qdd6GLGPRA3ciMuVAF4Vo3Z9MCZC6G7rCcAmsxR7LzIaUwIJizn/6p3/613/9V5En9uhKx0RwbJC6LZNKvS3lugCpnJLyVwaQ0idO6jdR5oc44Zmapkko0k17hWuxO+8iiz85w5ve9CZvta6szaQmkCaQJpAmkCaQJpAOFZBabase5FArFy9iiQk3aaHRvF/CA1XeXGzTGHSV2TDGoFmxoAXHhrkciD2eeOKJH//4xz/96U/5ZTJohfwlFzaCT5sDSBmeSAmsiVIYiYSfrOWajqC3LFgMDM/jprQDG9l4ZJeS3Z/80lMjeK0zJ+bLdRACURzpSP3www9zoCeffJLnqSBs4MeCnyWQWsJEdeqQFHF0hD08UgyDMcSpqyDTL50gXnaFhlpB3WKbZYBUwJnjlnoLxvCDMg690hT8yFTQLRIQY+QV+JUKPDbQDEhk577MjVU+8b5NKS0DsNUpifbYsWPoYzIU/fu///u//Mu/MAEKlIzt5Ur75r41Vw4EARn6hz/8oTCDjVrH5T6LYJJ+SnqgkDIPMtbYpXGZXMv0jz/++K9+9aunnnoKEbK/6EUvEiISn2bw3Eo5TeTiLfqiB+bwenuJ1lOMvQ+wqi084JA3ciF0MM9SzMof1KEiHKLJIrfccsvNN9/MiPVWIwl4AmkCaQJpAmkCaQLpsAGJmggjFgIPnhqQWozkInS1HU9Dar9j6TVsLHeoWWWsat59wqgWKeNRQdHTTz/9zDPP+HXNQppGhITv2JTVVjLxScWkYkK2IXxLklpV/rlNoTjK6rM0DBjjAk/bTZaWKkTxlCJKWq0HoTvXIGS0TVOlvfpTZXRKm0A5/Deg9yeFnn322S9+8YuPHz/eCqgBJE8Jm38QBB1xlAG96IiWWrje14Viyz7Jpy4SLTe+GWt7qLG9rV0UzPBUshiye6vsTHpofyKaQYT3czLsgVwe33rpJZCKkfDjOui28xEifRWApVI9GF2X94//+I8CJG5HNDQ1oc5ycVfZn3irSyKmt5599lm9j4a4NXt5hZgloBRTFUtnzTYZH0HXWCnH9AIqWKI0Lv47v/M7v/mbv0n/uTue6z7UbPPQJC1+a0a+lefLzwkK/WNMVy4Q0gcRSsdRBjC95V0wiUIbYIFZE/f48WgCaQJpAmkCaQJpAulwAqkYie1LNQgh5aq2Zua+TXH/3kXxqGXSpSBsA6kJx6atWxhsGMoXBRv45pfwUEZj2+i0IXib+GwfzRuQmoPGpEY7b6Okh0QFHqPnSy+9VABDxSxaniIjsXrLh8Y6n3qQegFu6mk7VBargIpX3GcYdUrnbV0NV4ANVjz//PPBL2ca6RcYKzcCWXoHJFwJVzyCIoAMRa5Zoo22l6WwbanP0n87rURDZSpUh8c3iG8Vt1fonB5uuukm4c3Ro0dHGuvywOMReJCFXbCqlXRez9JMdF8IuJpQmfX1Xzo+tithl5/RmIux/Xd2bxW6RwiSUX19JWcF8rIlwBtvfoUibESBVL1zY9DxiQI/JWHwNwJecMEFv/Ebv/Fbv/VbF198MRup1seStjdqr0m/eGvTqObcW/+2xHwbRVEgizCZ+tdccw3iF110kT6Uwsv4QVnlnDlIl7br0QTSBNIE0gTSBNIE0iEE0shsYIzmtYnt4p5N6c64qcnuc4jWsWTU1WnH4+Nx+/1pshTSVjq5YLxy/voQTnKGoVlux/CtWVodfNv09/3338+i7MreDMNp2gza6+Kiq6666vLLL+flLOQ+CkTrgzpQsZZApWBgAKmFT7jtsBDVmjofE8rLk/B4FcZ4A6xCEb3jHMHBZ71S068lK7Sdd4mq+O+jAvHbIbHtONtLfWwOPhbkDCC1Yt+v+5zpGyeKzkjkoFuhXvovyxZuqUJ01IQ+Bxq92zisOlYxqcWWwRdHlcEcnlsVVmpImwuUWtrkO2MJKTnuNpA4t6ZJzbglHLjQN1EXHz3nnHPOOussv0IR/R2XLTRFZ5wxswRSCanlxLSoDJEzzjjj3HPPbYlUubYuuIGYlsL9ui6kHAfv9eVguW6qzYbq+6gFGxCOnwK80l4LvUqlLbcVCpbZPBNIE0gTSBNIE0gTSIcQSI0FGaP1SMU/oWWsRBp7CbUkiWbHxt9jQL+9i1B4GBPNiotG5Bwd65wyw9MyMdoEpzzXFZCKtfAmyvJLNlZEwSsG3F7nOrTAxpSIQtPuTYkG1+ZYx97fcRVCEG8wPRjGKmW1YfdYGcUY2BZ4iMS0uAwSxlZHHZHiou28OVwT0MjmEOPInCA0tuP0O8w2TsoZMVI5It4ie8uinnrqqSeffFIQwjpQShtQpCcq4VLTugPvloQ6Vm4PVguT2ldR91FGRTkf45zCsUN3Z2+3DVBY5axUATZLII3teBj02LFjdYvqM9Bll10GAMDjV2xz5MgRCqQZTZf5Os72Xm2k3p/1etigByoV/tVX+pMJdJFtkVnQpcX6hdI13rMpw+7RLJuk3dvz1b58NLvdMi2uzihtf98K//6dMNbyk88E0gTSBNIE0gTSBNKhAlIfjxuD9kl+YCbAtCnkZxcF3ZL2vLXXYcxt1EKAd58onXZMwpIM+u5Om5QiqnFBQR61d+TyNOIxpm8PnY4GodYmtTsnpqyFtr32ShOmVNDmlU1flrY4pqozXjkBHVMHVy440IgK2v7SffW9C4pazGAlg45VQ8vdNlvU7SneyFgHIV4aqQljOXTn6pWZ4cKfmGn2ecSc41uChjwtffbJTXnkkUdanaUa2OhWKMFvHub17ZVdy8P/amJMWHemTovHWs3VgT2dfFN2BZd97aZ0OOJSnxl99DXiVWY9vilAVczW/uyaYynmG3tEthnQWC6+Oidm5LN2wh81ap0jtdTcdemkOuWx/3sIL1YPTsv4sFNe2jOoTJc2MAIhaOlLT58o+pIRLpodaEX60kYTSBNIE0gTSBNIE0iHCkhLlxpD9tAyThFels5gXkq+EnuYaizoZQAcUx8DcCmRTMo9uikuuBotQ1EqGGIvgTT2qk681hF1Oq8S/DoUhCAty1Eay5YcOT7tDzcqTGrtCokSv/NsRhzYeR6tVc6EuUKLfFZHOw/v1JYK7TvZBuXiigQcfDb7XAaksrNvGt8SSrj0VE1dSbsvsSttuNm+P3+0Kfl3SQlNxC9X+KxOd25XxEwjnhku2Lyz37BRJ9hCI8VFi6OW+hzHH+KBW3coY8ftBNFWIo2VQi2JLwjxZwcxbqNohHa9O/IMOig6lvKBsQNp23Mvj41cxdvtdoRC9fHctp7sTqVtaJWHd+zN6O/8jph5sDqBNIE0gTSBNIE0gXTYgFSY1MfsTmJWQktf3+9clJYil225c0S7TF1th+WA1GnHNCsWouIOHnbBw4xuW+yxjaJl+NEuM2PVTWbrbLYOKm5kn53KEkjmkaS4Dc4EH/qiwU4lLPlgHO7Sh3Cqv21T2hN8zPxui9+Eu7E7YQuTiFkUV3po/UKnMtc97eybRsAZYvsMr3LdhHfHdopjjU2qqIPIUsvvEyssqd9O3/TP3ZmgDNd0lQWLb7ENDH7b6zO1r/Q5khvqN1vOVPYA6QCm/T37ppLL8k7iN+Of4Dvj7b5/NAOuwti5u3OX2+rRnVJLl8e3fPBEWR3DrE6e2beW+p2+hQyf6Ty/ABwuerqy0QTSBNIE0gTSBNIE0mED0tKl2qFP6ThkpXnhUfrIvf0RehtIY7sWHLOB8eg4t+81J4qwoclQIjXrvULRUviRwqg0+dhXapz0kT4OlyIMmVexwfLD+RJLdSJljma2JB0H0Y3P5MNaO8Uv07S5/ps2paVWfUovQbP53MKzmtuOOVcnB/cZIGEr2aINuMv47PzF9BD/KzMt9TDWsb91U0qNTcBgufyG0T5EI+RYmWm551HA84rmcNIXhRWKXAtHB4pici8+y7fI7vVTGWKssxopuTu5Wtl98NmStuU5f+N7w9h0qRQQUqw8apCaQJpAmkCaQJpAmkA6hEAaHonRNmf80B6lL/d7QWinB7SrS5mXTNVZL8zmt+nacg6WY+6dQApLY8nKkufYHkwOEWJ1e9p3JXhGyvnqQbLZStgPLspeKFraibwdWEhYkvozSdsEMwCP7inE7vUtYZ+OIA386Ymy0sA+KFqdXafkVcM13/8fy3Yft+Rndc5x0doIa3VPn/zkJ0sAFYK6GNP9O2e9t4E0VqFHf/sbye0HKKsY/r27yntOlBF+5yHDo1b6nECaQJpAmkCaQJpAOoRAyipjre/SWUdZ3dkHPzsNlmHG6HYMu3fGRXsBabnT9DbPS1aXd3YaZq9+ZBlv7JR3LyLbs/91H83MDgcdq2JGtzVa3N+ZlpO5A8xLJewv/thtc5tmPd3OeGOfvmMvB11GdKOTGuF3k/I74439gZQs+yQr7NVd7gX4fHL0FyP35X0nyujZVypd8TmBNIE0gTSBNIE0gXQ4gfS8lH3cazktvtNCe410n0fe9qd50m7ilBoafcf2fPHOLuCkPK/gdBr6WWl4xer+08f7C7uXH6864nqrZeS5j6p3tr5PssI+Jt7HwT64b1lJsZPPCaQJpAmkCaQJpAmk/8JAOlWH3ic54OB0TrXpX3evsbLTc+kXbj+tcpDw9VTpnyqKTlvSpWL3h8oBUXTaIh9QhAmkCaQJpAmkCaQJpEMOpNOmeNJJ8OdL1/sP2Q/u+s/Rh55fkJ8UnM/RXqfK0gHzA5477Pdn/vbntZwqpA9ScwJpAmkCaQJpAmkC6VABaZ+5y5XWti26GtGenj0OGMkc0LN3TnTuEx09l9jgVEU4CJ536nOn7+7/zWA7BeFUQ8RTRdFJbb2PXU7PhU619V9TtD+BNIE0gTSBNIE0gXR4gHSqCzn+/0+U7wTJPothnvtk7qnKvk9+5AHh+uuYAT9IvPQ8zr+fxheOg5j+9JR2QEQ9Xz4/gTSBNIE0gTSBNIF0qID0fM0b/vqSBk7DdXZmxB58pvWU0hEOMiI/bfX+WmU/bXc8vV7pOYYiz0vvc3oW2Seb+f/FSBNIE0gTSBNIE0gTSIcASHfsu+JluSTjjhPrsVdr0fcC0j7zrTvZ3X9B+NI79+JzLBtuQ/DlBoJ7LWXfaYCxS/jOfSpPY+X5Hf9xFc0HTpTV1jx3/MdV9Pvbe5u9QTzx9xf8gI64/0KdU5J6W42tQl+u7d9ns9GlmDsXmR+86ZM6/GrZ2GprrT13EZpAmkCaQJpAmkCaQDoEQNpJ6wOLstzyZnsLxe0DSE4aciwXYC93aRl7xBwwcthGTsctv+td73rHO95x2223jT0ol7uK70V59BF3bHaBHC4+ELXU6V77Tm7zOTb3uf3E3jchvK2FKp1EcvtmF8KxBnsfIK0MtFJCu3B2SvEQ/KQh385dKVeWWu4rNDSzvzePswaX/UVqbLP1O++88yMf+chHP/pRF8vjjfeBULtQje1Et13/9hM7VNab7NxvaJvyUODYP8hbd2w2HO+AnLvuuqsd4VfH5EwgTSBNIE0gTSBNIB0eIG1zub3X8x2bc0TGoWWf+MQnOqH54x//OOERTfi9khqH8JXcnfeMnabdGT40HG4AeB8gDZ0Gns7lfctb3tLRvG984xtf//rXv25TxhnDS5dauVEydrjI8lSP4DcU2v6GfvnBXoKPMvZG9HvH5njmjrcZ5wf7Hae8qJxvtW3i0k7bJh/sLber7nxixG+99dabb7751a9+dQcEor86MXo7Vlnu6TMiwO0ub0Au5W8fRL1C5thVfNlf8EKew4U++9nPdtDLl7/85S9+8Yuf/vSn3fd0Z7Q8jgjK1mMX9dWem/GmrXaZp21NdzrLTiAFQhTSnvpvP1E6uRBNFuH5bVzeKYntVz5Oy55AmkCaQJpAmkCaQDqEQFqOPscY8fYTx3IYICLxuc99jsz33HPPV77ylS984Qsk10CH466OOh4lFjvHtwNd+DTr8ngmd9PTWhkO12h1Kflqljz3VYeOwk8HAcKMXxDqGJUOGO5cwE56pvoRlS2t3pD9rzaFglCOYaqkkFRpHN/5x+PU5FS5/ZFg8JnqOorPn1rPy/k3R/eLbQ2RlOz5FpX6db0TS5l8hFidiJyuMp+bpL7pppte9rKXHTt27KUvfSk4adSj7Q8MA+2d/9fJf0odyjhjZnRtt2/OeK5OIAkY23PWIdMj/HRMI66EbTh3nybh56tf/erXv/71e++91++3v/3t7373u9/61rc4GNlHpLQdHS0JUl1yjZMmFbyprBqp37wpHS9dH7pzcj8XxR6yXGicJcmXOBJS6iDbUc3+kTAQEbTVOY7/N2l1AmkCaQJpAmkCaQLpkAFpHNvWQbydUrb0MO8TkvDf/OY3Cez3a1/7GiwhiuJQ63J++fYTB/0CD+d+w6ZwI9fFBo10G4Y2w6ikiE5AGUZaDuibgHYHt3SKIOH5Je8hPF3Qmpr5ASIUpM4tt9xCp7DRaXkrZKpMRgihIG9hzCthjwhgQ3Z9h07EsF6dcchcPcX27Gp80gyasMcAKhOW7CN6wQ/lcFlj7s985jOCBCqlW76lFfxs+1NBQuccLyMrlfkf7XmFWnjAy1/+8quuuurqq6++4YYbeAP2mrVfTtkP2ZuGRsSF3+Z5m+rFnppFtpTp9SaC2SsXDDarg3ZqC7eaJiweWMd1yuQ2/Ad4SApRhKXYBx544LHHHnNTV4L4+Gaw6uWZhg5ZnPbcScOf3hQBDC9FH88YU1MdduwAbPyvou5l3EWfar7qVa+iN32QXzz7k7EoE50OXORaXilkqs9N8AmkCaQJpAmkCaQJpMMGpDFjmGcHoaH01IdX0dHXNwWi8ioU1cy6RSArcPJjYr/yla+88cYbuSbxMr/iohlM12GVz3FTNmg+txOLl0Aak9QeoQweLJRf0q87BOETTMVIuMUzP8PbmzalYXoKXQYJLMebC/mwhKDoAs90R3DCMvN3vvMdvo5Pr48QZXmo2yrw0C6aOIFAlCmdbvGAT1bxrjtEpkb6/MY3vsG3xAmc6eGHH8a8RyiE1VWA1GcDd/CmIwPyL24Kp+RMzByWuML1119PRV6BEB5WCDEcNH3mi3lhSQaFaozLHF4kIz0U5FBvEU6KTZMrIPmzEIuMrBOeqZT49FyP/KUvfQkRzRXPeOvuu+8mOPFZjZJHeLYUHMGAxOnJ6HWV69NRwy3idN4culf0m1QdDPoGsP0FpYBcTRpj8Ve84hV1x2CvIdf81n2aBLBrr72WICoTROt+Y3ICaQJpAmkCaQJpAulQASmPZx76UsmfKWIMNCmRGMRr+jumG4xipReVQLicryQtnjCHD8LjnsG+uikoYKJPyCojTpV8DnH+obligKWDjgG94hE5CclU2KPiEeqgwO+feOKJ4ZFe5wRUo626iZUzIch7sOQaJo8cOXLJJZcYImOMgz766KOPP/44d+dzoiYteqQa8xTg7QwSUpdB//33348It+vbObWgoy2avO+++x588EH4ASQ6ESBp68knn/QKbeBqJ5Aa7qPP13kk4i5yfS2izC7NgxeWEC3FwsASnEWSbnrLU4pKt3BVTOgXq1yCsKRmQQwI57CHc9XKfljNfTfpjzKdHz9+nCb1StwDV8UVI6REWRDCQemTB9IAPQMG11pmjYxEk1I6vMXTCEKBejf1sVTfynyFcIWLhf2lOOwEUj11IXcJE+hj/hOb4j62YQl7HPiaa665bFPgCuU+gUwgTSBNIE0gTSBNIB1OIOGS2DwDOe9Xg9PwVHKiwip4zVRGh2HAsJJ+6ULww2YrIDVIFSABkgoGlLT/w01hCfZDwesq8Azj2hyLNrWOmWW2ZVYvmGntCoZZF3t4LtgwOC6R1DX3feSRR/xyR3oJ0lrBUopbOlNAQpb2tWtAf+655xrZo0xSLv7QQw+xlqc86ejRo1R57NgxY2WVvdtCmj7kLwMkbHiXZ+CErzcn65GGCuGwB/BQRGT9i7ZgAFzh302ePYC0HNPnUqRgSAgZMaqbyd4EvZsqu4NO4ITS4q7hoH110HThLjW2eImk3mpuGgUdFgdgRF2qPwUkqBEKt6UaLJMGRnemIQoEIVjiG2RhINV4izt88S2bQo2USdt8CTIpRLsYGNPfo/vwWwZGoTVOSrUhlDp9BeG9w0OKtItgQ+AIuUd+DC1htdlzb1EpTGKAf7I4yqTDGO/ljf4lYPvyyy/nIaDR15q6kgmkCaQJpAmkCaQJpMMDpBHMEJtg2uszOe3wVEJeccUV1IEKGfgl7r1Muddee61Y4tJLL3XRvPaY/h4phm4S75ZbbsE0v/nFL37x7LPPchR8c2u2ocEbbriB/BzCoJwbkQQPr96Usi2HJxV4lCjgFShC3FPUMEaDRrSN47kFavzMW8Cc3qmpz/8taR6xR9O1farH6lWbQnZ/sn3hB9lRvvLKKy+66CJP4Z9oXmR7DCPrerkcnd9giVv84Ac/8BsnKmjFNSDdc889PBKEStP0Fkn9MifOPUK2VNcVkFrRVADJvViNa3J043ji0w/imqaivljQ9t/+7d/+/Oc/H7HHyGxAHJ8E9Iq24AQpvQMioIU3rLKFm3WFWuRbAPnYY4/5JQuH8dYy/BjfEgjlvnd5jjvE0RanPP/8888888yXv/zlId8vnfuTSuEkhWSdEca0jA3lj22K6zy7JAOm4Qy8hQa0SIGta/IW5VCXOgOZfdEpp7aYTWch0CJOdkS/DB7BpxdRwHMT903WwxI3cLPFYxllAmkCaQJpAmkCaQLpUAGpue8SSTWQp+KJpx7dFHbi017rwznqAicudd555wkY+mi92sxlxEgtufGiaOHpp58WBZU9QMUgighSKoQQWqCmRqLus9NIn22u1rsMT2CPrrvuOjbAMOviCv9YRRPIKa7JWTf9iQJfoSz+oXKTsNvbLMI2bstTbJmQwtJ99tabBHsQihOkaL8ky76RL4OEZuH5sUYpvUAie+QQ+TqjarfFWhpFxLsFPy37WboU5TBESbotasLPNddcw7R+SUozQAilfZ+AKAGSzuuJJ55gAryNqedgn00xAz86RIYmrzvciFx0S89aQdmLWGqSPZFxwvRMgPlkL+DEP7UUdHnkxXJWtAiQfEbnS1j6gW2kwqqCKy9qJc0vw0Law6pQEx02ZQjWx2om0FAxkvpELrzkpRxYQyV/lARcwsRIxSapOPZ73/seyt4lEROo4BHTaBdZdFzr4gWu6rtTP6tRXjeBNIE0gTSBNIE0gXTYgDTWY5RJyZUZHruu1b5hU17zmtdAWrsIeQUtGqd3enGdkH1LXi0PzkXcJAkIYVfzrfYxosUQnRoflyJAL8yG8oUXXnjBBRdot3fHlGWf4fuMTe+8h6ZaqU4Fx48fP/vss71rAN2KEQy0EIVyYfjHP/4xHogWq6ttAemhKdGxoQyu2ObGTeExlMDbMqcwJjpt+lMaazkTrdz2iNJVU5lDt2ypR62Sai+h9jzSUMvacesOUq3oXm1RhII6rTUibHGRwqWgCM6Rop97771Xr8HR1ecueBAD5K/4J1pfFIITTWIP22RkDqho7rtIGEF3sIefT25KORMldJbu4KLF2COlg+bLxG2+u7XlWsEP0xTLaeKHP/zhU089RT/s3lcQr+v+NNGHiiWQys8QZrM+Q+h/vYI+DDz44IPghNUlkLDBQ6hFW5hp6n/sZ1p+TNFUOR8swt9oFZH0rxpQlJABbGJdnZHW6f/qq6+mFkQoeQJpAmkCaQJpAmkC6XACqR0bPdOqsaAhMlbKWG1tBsGMv8nZnKzS7GHLpOGEnGN5xpgFbrc+NRlmLA0qNbbsT4FQQCWAP/mEIOeSSy6hAk23meBIGmgKm0bQ9OLFF19cHNWSp2PHjtGvF8nfrHTDYvWFZ//wD//wi1/8gpO17VE0x1dzGsCkVwiuvrZ4D67wJmCjAYou96IsSTYoK5QgLRDSqIuxKWR78bQ9eurqTumzYHl8UxDPEZscpxPeqZfJNUf8tlwRzRf1X9xRtXKF65Xc9DomywZWTSuethKJa7J9bRV2FtTl+gVIimrMRA8qqI89MlJjM/UwVk+nLTopeZTbjbnvxOQSJeCyQm2VSErt6Uf4xCl/8pOfUCad1FC5At5tu6X6teGfrkHu+9//Pr8nr16DRe6///6f/vSnolCUscQZSEd89DWET10DUnpPeCuXN4S02ymWWl2PlNd5HYvQWDBDoQxjyHzkkUeefPJJ3DJotuMebVqUdSaQJpAmkCaQJpAmkA4VkAoS/E1yKPKOwahrMmc8inbREu7mAR/ZFJJwzSKBAaSBoj/flOYZPeWguE8G13xRi31NLw65aVNaN06eJmeXWQj4ab8bLitOuPzyyxkAb5TVVkFtda0arfGGNMX7f/SjH/3yl7987LHHsFHG6jIhsuQJj6ib7ATEZD4NS3qTPKmFQzSgMsfVOoVq9B2bMpboB5uWZ5POLw2jQF7aaN5fdFdnwd1bG8+BaKy1W6UXjITLVSzXKnd3UkUr2LGKJc31MaPJep7kUR0fCFFs+xYVI9WDtPeTdnkh5atZIIQx7KHpmqfy1wceeECLKJTVwYiU3EaWY3OigkA3794U1FrjzQeIz6MQF5DAGM/RZ6lfv4A9LwIAPbTrQYkXwz89Faj8/Oc/xwYlYEMr7PXQQw/xJeKfd955LNKHFnSIdvToUcEMrryiOa238VNob4P1km/wQFd92Oibgd4WD1rRVlPkCrKY4RKaazukYqQJpAmkCaQJpAmkCaRDBaSM1PLydnJ56qmnNOz9vscbFrcfn+bvuecervbss8/+9V//tWpGn+1c03bVfYkfY+WGoUU1/I8v4sZAk/B82gVLYKgzTsbMLE5KfGxHG8YbQBq7PBPPwJpIjEQjJQq0OGSc5dbMewmsj22KQA5jY4PI5ZJ4TkyQp59+mlwAg+dSPVBAluGNzt0nKV0zz5lnnnnuueem+lbFt5In8RGEH/6hJl/0p8AD/ZHHyctFgy3vKTCjQKySiEVpu40sV0BCpw8AzKT/onxdGPa0hRMRDp00C+8OrRIWw5rui4WnuC3BYuwHWmmljXaLRrzeDLsLTfzN3/zN3/3d31EgOzYrXT4uJ2HWdhxYJoO2akgT8Ewcv3UlXveiPrc8XXprq0dgU99N1Ji7I7THdgDJroLQlHWEu22T1IovbYHBZZddduGFFxIhO8IMeQXM9Iyg7o+6tFiANPYYRbYNRumqTZ2wx0VR8IhKQRFmUiO2S6EoBzr/xENeNIE0gTSBNIE0gTSBdKiAlO8yCWH4Nx21GAN1zdCpMT078XUsCp8effRRdYhdzmUpgGMv7GXiZgmabEkSb4lVfvaznxlwQyM/wF87fXulXAT01ecrfl27M1Zx5/FFSq2N1wrK9MuinQiCGthTKA1y4vGxn3aa0DeGhpCxOdFyT0PM6ESeeeYZGsjjycvbWrhuOE5eDADGkSNHStXFQzPp4xyaPsZzRy5CaS95yUuuuOIKpNremkSedg6Nm7gtZiAF3rDavHabIi23FE92xDFGEPJya/1RJi8pQZDTZHfriyi8deyt7ypaW3YiGH73priJpVs2pQXVeKA6vwTUhGBYc/rBogIENa3ToZN6hwGkcUofglxQ0HL++ecjxV/LBO2gR0Q6OKd5f1bGMLlYqk2pVjs94ZmedeL6bt0ZOnHudXJpCOeko8NPbQq2KeTss89mI83xEAhpjf2I4ds2iJnYAj9plboAhkp1Ii9+8YuZj/NjhsiICJ51/W1rzvP7uDKBNIE0gTSBNIE0gXTYgDQ2Cmx5RkNY7ktOQ3YUL774YvHAJZdcwswGiH335T2dP9cxeGMTmQGktiJqXQcDsC5fbwa5TXM6MqQhdVvMtCtfCY53bor7GWk5BdxeMHkhmsyMLDO3/hyTHB3D55133rFjx1gUTUJ52k5Jq8Pbmg/NVPQlJCgS6LDkYrBG9qVgFtERGW+drteOS8ssBDY2cOdGL3zhC5lTGICIyiiw2RipF++1ySb833TTTSCh5wrqy+MPI96ORZghLKP2/R5xNoIi3GoFfjj905ui54JeAuK5jOSlmcb0d0kVJelqHRsv3xTM8LO+7tc7kBo14HloU3h2Ozcl0YgPeY64CEtMwB3pwZ8dwTJOQxl7o5Mil2N077Yoa0Sbg9vMzXO03m5HbbzuaVuLdpxzq56EN2eddZYYieDUnvXH6XpDfM7MlORq1dYTTzxBsZ56Hf7POeccnSA74lO7fEwPi9v2vsxwYy+qCaQJpAmkCaQJpAmkQwWk5QoNY8oW2SKBIl8UD0DRpZde6lr4YTxK3ePkY6wM11xu0j3OX+kzf0tKlLGaheHbalmFNmRRyoT4s0UZe0APINUWhrHK79vHkDbzfuY3KL/mmmuuv/56DRG7M+GQbRvowo/BZ4P7vCT3pV8yoqMfKVu37bA79aQV5sb0lKt1PcL2Dt0clHT82yD7zDPPxA+aBEm6Fhq1gr0FSFR98803j+MGV6ceLjey9HqLoFCjSXx2sBwinuoIRm/FJ9pwp0NQive2TzRxUepJa7T6nN+2QZDQ0n1FwClIZnFvaQWMC/naSXx5xnP7xXeoDAjxZnoo/7jDjEs3aXfLFurjszNg6pfrR1ZACngJ3spzf1I+MYWCRBYGhyJq5KswzAH82XE15cuORfuFiG327RqMxXtiJCzxkGWPhjEUOi6auTt5sf53efbhBNIE0gTSBNIE0gTS4QHS0u9bMe6iVd9lD7z61a/uxLIWopRd6n7HG2Ml11wCadDUWF6FCMOwB6tceeWVlItOGww1JztUNrAdr+PEvqaA28Q5VXaGR9Ov7ZzYAR5Kh+Z61OFwVKlyB6IMzA8+2z+o3capD6vXXXcdYcnYTHTLp8tOaBfLTolDfByUsmS++V9aIqwIk7Cg0kKguoN6E9722k3hoOP7+rYqlt1TwVi7ebYMHpOu0Rz7E/HRoogSVbEx8Lk60WR89uCCn9+UMpXb1rO0YNGCgLOFXqTAgMqCimaxl4fZLNfD1ztoPY9MhxmrJUBNHxdstwb+Dzcld9o+f2VsoqSC3qflYaVQj57UUyoFeGr3W8JEyTF9QVkCqf2h2tqy7RXaYar1bwXtcFUKrIb8Umz/FVZnNE4gTSBNIE0gTSBNIB0eIC0PmQtLrY4eYUOrX8ahtrduCgu1mySVdd7GcpvFJU1/ep1OuWaLcI4cOWLkbQANh9RdCmDAeP+J0szvTlW28DjMlCda06l+bPxHNVRQHmRrzjv7OeIrwLsmBQtl/obIRSxFcc07522ANA6ipuuWSy2DhHGyNRUh1foi133IV1zADB2+flM6m7nAYHn43+p050oBWFsA9DqeO1ilLqZZbLoVOBU+qcNSHUG99KTlEi+6uueee4QZfnkPedmlHX/ASQ+oJ/VnqaKlnoyMk1XJQbGBwm2bQhV9X1HGYrAmypOldIoSl4dNl4Iv97pqt3EXdQHs0kEv3g269fXqoJM/Lze3WvpSuR3tMtBxOMVszcgr7XTfPublp6764gmkCaQJpAmkCaQJpEMLpIGlNrtJ4JZt9N16eDAdtfymBTPLD/A7Da8mHbVLNTgxiaiA4jq82dOxgmiJom1VLuUPSDgZOwkuIdS8anppoVSz8KvdyZfEo9neRm/eFO7YxkClFPRVnrpLi3SBbEYqQFrxWSBERnQa97cDZr1S95U6oyaml9reCaQRJ5RQi9X2YwqfMU+xUFRfQOdNfI/l5TvjLlIIeJrQB6TR73QaXwvpC/DaPoD7Jvhq/8rlnu8Jm2XLRKlyj/CTI3VC9tgiMxSterql4F4ps6Ss39zai31QecOm0EAie+tPN2Wpz6Vic+k+wLT9Qcu3PrMp7UTflMHw/+GrE0gTSBNIE0gTSBNIhx9IA0uZs4tw1ZzgyDBYuuNO86x8dMifClJZPDW4H+BZUd520CE/TRWejc1rsvFdd93VcqaSX0t7LWFxLz7HVjUIZmBlDNZz8eiD08jAGDrZjuVWW0+OXcXzhtxopF3mZ0Ptq0SBbSUMVqMzVuYUk/DI0m23g9i9gFS+cuejKEDlzyFgHyHyhw+fKKP72OkA40NFZbjvEkVtSNphNgF1m88VkNJn6Qh94SjgdE1wwpbwmlZD40p12yH32EOqza2anS93JA2Us5LOl766EnkCaQJpAmkCaQJpAulwAml8Ph9lKDEDr/xmr2/w2z46rL6MrE6KyeU5xNt2SlmDSdalDpHS2Jgyky/dfS8gLQ+Qbtw8jDpyW2ti6e7b2tjmc+Bk9EQfOFHGUSgR30lzr94kJxhlCa2B0qHqnSIvl9mX1DLWgA3VjTIQNZhcBXJ7YWnpdt0ZGn7PiVJ/uk+wvYq+6kPfe6Isiaxca9uXVh419i/o289wmOUWoquyU94JpAmkCaQJpAmkCaRDC6S9hpX7l32AtK3WvcbrB6S/LH+8Vf5kUZYmPymQhtaGl28namwTXPnlXkwuE3BXfrwP2b1aWSp2VQbz23rei9Q4imYEw0vVHZCx7ZyJwc9qO87B5LLs5aA7+6aVjNtevr+PbZMd3dleHcRJXXECaQJpAmkCaQJpAumQAGmvcfNJ57X3sfFOXsfpv6cKodPj54B87qS81xT86fG5RMsdz2u5/cBlp8dvx8YHCfxO6gbPSzlgD3VKBFey7y/RAU0/qk0gTSBNIE0gTSBNIB0OIP0fj8t02BByAlUAAAAASUVORK5CYII=",
      "text/plain": [
       "BufferedImage@e997703: type = 10 ColorModel: #pixelBits = 8 numComponents = 1 color space = java.awt.color.ICC_ColorSpace@7d546d71 transparency = 1 has alpha = false isAlphaPre = false ByteInterleavedRaster: width = 840 height = 168 #numDataElements 1 dataOff[0] = 0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val onerow = topvecs.view(28,28*opts.dim);\n",
    "val nc = onerow.ncols;\n",
    "val tworows = onerow(?,0->(nc/2)) on onerow(?,(nc/2)->nc)\n",
    "show((tworows.t*500+128) ⊗ ones(3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAACoCAAAAAD/IJEYAAA5m0lEQVR42u3d17NtR3Xvcf875sEPBJeMEignlAMgoUCQCJJIxoADAq7hXhDRMkEJkDAgoXyUc3Ah4cKpbB7ssp/8l9xPrW+dvpO19l5nn3Pk643d/TBr7rm6R4/06z1Hz9Hdv/P5/7py0003HWWFUW2U/2yWNnvc0vVNu5cjYOamwy97pLkXobrfC+dHI+/oZQufa/zsJsJuxHfrbkc6e/ec35lAmkCaQJpAmkCaQPqfAaS9c7ObEvdus9ed8lHi9nC9fC8mPBoU7VH5r0sXe/fLIyC7CfIt/B+W4JtOuyMmJ5AmkCaQJpAmkCaQ/qcBaUf+9uIQh9XFkRHZ8g69mz2OMhI4LMDv9ta+RZ9HrMA9Rqo7Ou5eQpqbjq4cZYi4o033Prbukb1uJpAmkCaQJpAmkCaQ/hsD6ZAedgTv31vM+blVOZqX7COD5VGGJVuI7Gaw3R4eVix6NJHeEc9xb5Hlc4tyWFja0uqQcHrdo9wJpAmkCaQJpAmkCaTfaiB9YVEOK37YTYN/drAsVXPIL9ObDZdt8fbFRTlct1irsyny0aNoEPncb5a9A2k3u3zxN8ummXYjst2yRwyk3YhkwT/90z/9k1VxsxR/tx5HK+XPtpbtMFsTfOhqTfbdhuwJpAmkCaQJpAmkCaTfIiB9YadST//rN8umnQ4JJNyku89+9rOf+cxnPv3pT//RqrjxxHMq2DFJcdMM6PzxqmSG1FdlvH3pS1/636vixp/xuR3bazYYKtb8z1clqT1c0+whNbmjc2uoL8z/8cGS+JsG29TnknLsEfPLX/5yIv+fg8UTP62ZaTcUrXG4he1N2G/JAFjjcxRPVCNybsD6A0u7+dIw+mcOls8uymcWZelLh0RR3qJk4k1JB3TXOJxAmkCaQJpAmkCaQNrnQNoRMImdnTIVC7m6X/roIeeF0wIJweaTn/zkxz/+8RtvvPH666//8Ic//JGPfMTNxz72sU996lPUpOYh37AHKU1c0522ZMYP9r72ta9985vf/Na3vvX1r38dw/jcAiQ9Doce3qwCVZAdta985Stf/epXXUmdB6TcTbUuNTks9OVF6Un8JAURPrkqf/iHf0iQtWhhU59L8MQbYb/xjW9882D51sHi4c0336zOmpl25DaCw6XG2LF0fTXHULhF/PEwsqRmAgrEJ3PgyhXbiFNj/kB2qlgOo5vDXKijrk984hMfWxVe9ImNEjXg3BFLmyiin5tXBZNkb6AMsQOinz5Y0FzzzwmkCaQJpAmkCaQJpP0JpKU2l5hRciM9pYiK+xxraaQt35VxgA9yws+HPvShD3zgA9dcc82VV175nlW5+uqrPQEnaiJD/OwWdbhBjbKo8sZV+ehHP0qnGYMUWP3Od75z22233XHHHd///ve5FFlS0Ob0dyjC26dWxc0AJFUgpflfrAofpYfMH3SxuqMn4UFbqsvFR0lviOTc2iJCJzfccEMDCokwsCOWel+PqwGeb3/727fccgthv/vd77r+5ar05/dWxZ+q4YRlN1WaPkdI+dVVQXwtvsoxMjTeyN7gtZv4MawVIghiAJOYySK33347o9An3SJIXlKzIPMxxI7DHIUzik7T1Qc/+MHrrrvOlcb8mQ/csCrXr4o/1dxxVFozE6kZBTM0iVX8aFUviKMzsPrxVeEhmKztBNIE0gTSBNIE0gTSvgXSUpvBZrzO9s797YMll6IXV8/porBhe/oCedIXd3n/+99/1VVXXX755e985zsvu+yySy+91PVd73oXUMESXlUuDFgzfEyyqHv+R8hrr732ve99L0C+733v0xZx78d+xdidd97505/+9N5773VlOeJovpzaXoZbeCtmY04UqFI1fdEDSZkfNR7AOykEA3m/4gaFNT5DkbbUlffceuutXEdzbl3YRm/0rJr6etd1I4urewyszYOP6WZcaYsrBFEmmqt7T1gqq+lCR57/4Ac/wLlf/cSymWmIvxw6w3whVgGM53goQsi/E5NDNxrSWOLv+HFisIoTqvvxj3/8s5/97Oc///k999zzk5/85Ic//OFSmY0j3CMfXftWERtgxjqqURFz86KaFF5ixq8DSG78BHgD6rsFSMzEkw09BPcQEc55ySWXXHjhhTyTdyE1uigqw0xt43MCaQJpAmkCaQJpAmm/ASmFjjdvFHvzpoLxqq1892Bxzzacw01vluP9e8fkwmYtcUByKsDoFVdccfmi+FOY5DntkEST4rTmiEfxZ5GbvryqUij4XXzxxcSGSRRqrndc3X333Q888MBjjz124MABWPIkJ+ileWl17PGMwjZXfKpGIXyLKoiJFCf4q7/6K05JLRjg+l6gVaZf1JapA25oIwTyY11ry420pToEm2ktvOHxFE4uBA0K7373u2kDG/gZnwGW+nSDK1b/0Y9+hB+M8VF+gKXCtgI8rTCgGn/ltXfddRd76WgtTBqBXHbHD4uj5t5PSPGha1eFa/LjQkduRMkeEl+dzWg211d0hyYeDGcPPvjgQw89dP/999MG5oGcfiCWQekfQdS4aWPoAHyk9Kv3QmtxkSvGdK3aV1YlzHuSx49JcKwm5giPN4GkeZ8KVEP5/PPPP+GEE97ylrccc8wxp59+Op/UdeMCasHJTf91IjuBNIE0gTSBNIE0gbQPgYQVCm0ekL15T+/f7FR4MN7Fg5AnvX/DEp8YYdLyvXaZmYoDirjuuusASSTA492PnAYlv/SECtAZClo6aIEHL9QXpfO5c8899x3veMdFF13EBa+++mrWZRJisx/ffeSRR5599tnnnnsOopiNkbw09/l80MQqDjGguSANVyDq+YgZCMhfvdzzBjQpgZYQ0ZcmGmb7+BxfyjFAP/zm0Ucf5UNcmWLxn1rUR4Se0yre/EQ/xgWyGBFylPgc+lQ0xBJHvO+++wQb2uKTpemzsYnhuRodUoIuGOjeVeHNfCXbr+mzSLiwKqgTgRIIaHg677zzLrjgAsRppiEG/26g3RP6bBQeNEe6hproIwjwhUa0R17wbvhoqh3/DYtNWJO0vIqRGTyi61wl62CYBxojjFO0QYE8sHn5ui7xYom03fY/yvOJTOeUf9xxx/3e7/3eG97whje+8Y1nnXWWTlFuXgAdNBtQmk2YQJpAmkCaQJpAmkDah0AaticbhoCE9n+4UbgF1pvJbV717lXpHX0tK3SZn6rj3iwx2jsu5pYTxxVNykekAqyXAeBmGH4ELYp7RERHZ599NjMDAASW3ICCyvyY/TjxCy+88NJLLz388MOMl7WWmROln2Ky+IQHM5UnJQ2UDUB9xg4G4w3e74lMA9ryJH7Gg5sMHUAqGUId9Z9YFfb2hB7Irj7OWYKAvOrHq0LhulMBkAwNl156qZqDz6HPQkeuiRl+qQvPDUyanHrqqaeddpox5bLLLvOE20EdN1UnziGZBvpysKbPvmSM1A06xx49UMgZZ5xx/PHHn3jiieeccw4seYisCvptwprsZTwsB6ZiCfz7CbccnbfAD02qMxKCg3Gz7eU3KEgtM0dLOslw7KstOsR58sknn376aQGwe8SZW5NcrpSUUhbKqEioMSm/lqZc5Mn3+AAgUaMY6Q/+4A/e9ra3sa9qPLwpe0Mq9niLVgWW+c8E0gTSBNIE0gTSBNJ+A1LTtbihd2JnXaoPP5tYclXtJ6viuXfrUjmXQMJlKYbLFSN9GG7ikrP2Ul6q5cgr7XM4Zgq9BpA816TKGIZJQIIigQoXLBOiWIJOVQOkAwcOCJCef/556ka8j+K96I8RxI3miHBfV7w1oOChYKAn5BVpiEyEPaQmL0icf/75sBd04zP7YZ72wJixNUEqbsUbHJ3BUCYv1fFyfFIsZ+J8LAcMl1xyiToBaalPf1IIcyLOUiTl3Keccspb3/pW7n766afzAyKM7wdlFTTkAa0uMnz6pENyFQ+r1rx8mRalGgDSySef/KY3venNb34zoPoTh4yCoFZYYtBYWuITV/TZ6jI6KQODQtQpfaRc0sJg/RaBI1uKSdPxEUxwxmoBkj9VfvHFF//+7//+tddegyKcUOOIpdtuoNVouZCuCV4E1ceY5UBfmvJIzsBV0/okPW9VBIHUIsSFW9cGUErAkk6L1SeQJpAmkCaQJpAmkPYbkFo0Qsss1Idw1gUk93esCna5EWOwX1PeruUgKh5Sa0Aa349HlmHIwe5yKbjuyJZrNjGaM3mIv9aH0EKfvQeQqJjWmohXR6RRnivB9F7iBc9Iic1Zc2JYeuihhzTBwJWrwkXGB+/W26DgeS7uOTZ0hE+BCpy4qtZcNmpiHgQZ/pprrjnzzDPBj4DDO/tSjhNCPf7442Iz1AD4rLPO8s6t/uWXX04nZAkPBUgNDRjDA3AypCYj4XKsq++DOtMKZvQFM1Dkbf7YY4+FIvATHXG45tZHaixWR/CTD41YjpKLhzFQlFtRgckICD9vfOMbf//3f1+MhCXUHnjgAWHnI488gnipvWtAYmLuaMAiDlLl7BbDqE/V5MIGUryc+ByAy9EStok5JqOXS8XqSH3h7r/927/9+te/xgC5kGr5Vv8Myr6l/wI/9e9cFfcINie+li3RIEhXPARO3CBlFOBabEGZnrP4U089BRS8SxNyaVtycLPqE0gTSBNIE0gTSBNI+w1InC9aHJoHaMPGhUMeUjSBSetXzwHA1Z9qloAoIFmLkUqgLCjSd9+Yy0PFAeXqwjvus6uSayIeINWh5ZbWjC/xrp6z91+tChlQfu9730tm/WISPzhRoZRQ4pQ3AACk0BYeuBoFaRKQSoHVBW0CkvhEJOMhC7mqRpUXXngh76QvsBTz4FbEBR6spf5JJ50kLCmXYgQezX0bI1TDEqUL5I477ri3v/3tUMcMzc7nwXgr4YA3eCO/YlUK5MaWOiPgTJPpBHvIogmfUIqZcjIwADYt5FYz+CVsgCxbYuiTZYmGn8bBsluIoDKNnX322TiHUkOMh9z3H//xH//lX/4FltCPw7EjUtZnd5ADQmphQfSpgsOU7Ju/hSJ2NNI1f83uASajrMXbShj+53/+53//938Xeer0qquuEnYa6bgZqcfCOYyN9fZ64Wwk6uNKQCpA8mcr28mODTSNkiyCsjGaPrkKAVmQ1FwdZTYqXzYL6q6QcgJpAmkCaQJpAmkCaV8BqW/JLZ4BjNZmZ+OSU1uNE2xaA+PaEzVHbLMjkHII/ZUOSzb1+zL9yiuv/OIXv3AlAyDpkS6aT2zqdrwrZ3hmoBcyUyiQoE9g7oUyaYN0q0rGHs3jMzmVcVCeJ/zggk1/l8+hAg4ZnuL4YrrGufrnnXceIPEnDkEt+HzxxRefe+45faEgbDjhhBPYkowj8GhOmd4ad/DGBTlikGMDFagrFKlWhigGECltd7nyZxlwNiqRtz/7DN9cuQGF12rCBWmgnYN0hOb73//+cliLECI7gNRHAl6OjaKm0dyvuuNPaYA52AiK/uM//uNf//VfjSkqlJ6yPDiHHnBIjeI3qCagQJHG1IcoVmuBfQm7hbsl7GpFhFJg/TkIDq1iCal/+Id/4DCaG5Ig/MQTT0z/Lb5Cs2iQCdrD9GurgmA7ELku0264UEMkszafjgjli2Mplj+oBgX5lbYZCJww0Jw+JUwgTSBNIE0gTSBNIO1DIPWNmWlhg5e3xqPFSE2Ij3nq/KPl0yODc22HmiWQWoldxyCked+JvWTTCAHcA0af9rUq2VGTsShlaJNqSCXw6JO2h+Uo5Lh+gkOSt5odMFo/Q1+92vIJSvdiPbaqCWZsoLtWrfdmr4I/BUh8VH2CsND999//1FNPiZG8MbM9P4YN8QlkjqXII0ZSv4XlJBJ95e4QQu/4pDReRbdNCo8l4q0Sb1/FpTL7Bo+xG2+8sblj3REND4HEk/Z2LNe2CEd3HIIU5TroYmxrPoBEkIJhDjAWsbM7zt3ri950qiajv/TSS//0T//061//mu1UoyW/ZqzlIShtAaBTPfIilV999VV64zxatWdQzLTpOd5IQRZuWsrL0u79iRks8RwmwB77vuUtb/nd3/3dY445hoB5dskZKqeu9v6Ot1ahE2c5/d1iOXrjzOLeZ555hskwxk9OO+20U045hWUxOZKksZq/UW/JPRNIE0gTSBNIE0gTSPsWSMVIMDNipDJWm/huX79W4wSk0lXHV/mSHZe7CI3Mht56W3eNLF14XX7iiSe873rnRhCF5l5BKIUKJEpqHASXQOLHYhU8FPnkWH1g5lKQwGXPP/9879ChosRElJtZhq4mkccI0iY+3K5tjNqkRmgh8Gi9jQot2MYwzNMgG+jFC/oZZ5yhJm7HuuWxK01rtAjFpfTLSzBZ7KTQxsj0bVFWC58GEpZL91shH8Ldt91P3wnaqwjZwg90MKwarjiEmJCfJfJYfb2cqddXG6M3O0ylY6vElve0m5KHbH3gwIFyN3FOrrZ8ovOxbj+7B6Q23ylGYi/awzNqjVCiO6GXq6HNn/QTKWyUoLAJpFwo5Z955plQ9IY3vOHkk0/WUKsSYvzUPptsx5Hctzl76+LGTuXLZBH0xYRAzsToGHqgCFkuxPrliFBsLLnSTB+EiEZ7O09/TyBNIE0gTSBNIE0g/ZcCaSyhweXIbGg9UjPdLWxu7rsP4f1aGNCS7LXTjltG0hHLrdxmME1GZOXP9n1uc0bG5scloZavuTxobUx/k7kQC0t016aBCjelUFpohVIfvJvmbuUzcKqw3FZ7jCClxjI/HtoJkcs2X5xaWtMy9v7WSrVzzjmHxjkBykVc4yy9/BtlzzEGz+r3Hp8NmvMtv7aj7PzUqSf66hSWZcwZkMr9VR8n96xKdjHAsQtfVw3zuvNy39bV2CNyA8c4anoAqaXmbYkOeyV8tDlOXx3GdvClgQJY2QmUw0x6QRyilptrI16KRsm+KLTjpJ88hByxx/HHH3/ssceecMIJ7s8991xGb0PMAuNSYNd222yRUp6Gz1NPPfVNb3oTOgRscFEaiPXblgG6Q7Oxozh5uXtU66Y665EeggRTGn0aglErP7h1cSUBs37/WhoHA9gE0gTSBNIE0gTSBNI+BNJYj9SO0sU/LUAKTt37daw5b5K0hb7LT9GbS3l7z84YHQPThoA00jwpybk+wzO/P4uRloeajPUzTMJ7Hn74YUyig4h4xguxl2z2oIV2SOS4ZRuUuOCGXjrXbXnYcyNI5w42b9vCYz9ppa9m/wVm7jvRWTVEWoqMZ/w3/7vjwTM8WAU+VwCgu3a3DGbLHbdpvnVftMpCYXsAqdnb4hw/0YCh5Lnnnnv22WfdiNxAq/0lW1V16arkyq3VGUvF1lTqeRO773vf+7SlzE6DTkUjXAmEDY7lCqfnMkuWxycXa6HTbpUEZ2XmwJV4sg8GrsYgIxEm9dspPgVCHRO4PM5nMBmk3RjmNGy/JAK22TfjFpjplCB6xGdStPSoEHFYKqEYBXstmO8AyEK+VuCX4dFJfp2R0y4MbaTeYNeIM4E0gTSBNIE0gTSBtK+ANHYR0ubWW29tYfnYM6jNqUe5dVXqo5OJdzuMuRnwvLNgKW8oXKFo6vDSLKRpxxZAogWvtu3SsgakQhpu1I7PeOOR5a2iBn6FIulF87Yiamcc2h9z8WORzxC84wmbvu/P9vvWRXEgST30K7LteY15WNKvP5eHBy/3v+6oEox1GCH2OvitSdgSC/Rbkmi7Mumxyei1Y3JyhZZbsx/xn3jiieeffx6QHn30UbhiOD8VIHGjjrWGq7aozjvXdmYfPGNGTahjDm25ZqmxfVSgRpSbPu7Tft8GFDfj7Jnl8cnNQSe4YnQzPmbiuujI7QKYRnAij0zTVjctz4kZ+xc0Fa5VQ3DHwChtPMr6bSjPJXKhviuMvTWXgC9G6iyftlkvVbq9DPzbYA52afNKD0v46L9IK+eXG9NPIE0gTSBNIE0gTSDtKyCNMIkkTKL9QEu7MbY39C0HS2cw62ATRWsn9jUDTqpmVwlZPEOhYHP+qnRqiHijJdZq5mrLbZqX07WhvcXt6nRAbyeCUMfw0VDE7coeyBfHp/0lOJuALmJp0RTBS9ht58RSHzRsB2quhs8cLuf+/EZBX3f5UzPynImztrcOnTTV26acwNPZObrrtLzlUYUj4AyxmjDkPffc0+l3YkXNy1ItAVSnIaH9fTo6ZUcgDRvFZ98e2Kh9eVrNVeJFi8ADUt8bOrK6vIGljUZyAweFakELygKh1IUOqQsLWbAjrjsZqOEjFC2PPxyhV0BKM20A1Ml8HZDcZlVKuaplroyD+hqU1w7h7oAWbkmc9hwfqdUjX6fjIZe4cHUfq2NeYAJpAmkCaQJpAmkCab8BaZzYN05i7pzgvum2XvfrB0tnrY2ZyjWxl0bqo3IGaBqUWgVC3pvHLG1z1p6n67SwRNEQfsRycdVJby0syZy9Q+OqfZnx6doc5Tjxbi3uGh/OVdOqvRTbg8ZNOCx8aja/QUEZ0c5agDQMPwKq8gxIWmpmZ87hpPRWXeir3S3pvHSHzZgz5nuP12Qc9bc8DLsciFYcjRgsSy1zBdaAVGYJ5Rvjil6aZ18ui0rJLZr6o1VpinxTn8newvgB6XYaYjLgMWSIPEuXKYGmDRzLe86XNt1pACmkjbCtyC1+PFlm5Q46Yx+iTcC3F7miucqdMbO2M+n4RtLHm5FSvLTRBNIE0gTSBNIE0gTSfgPSWPiRP2UANwUP5QeM+1ItmxvdfJ1dnnS7PLcvzXp1Zi0vzczmHdq1bW6WU8NrKFoaKVXGw3if7rN3eKj0ahu34112C82at1q7waI9jEYX+cc4W3psYr5mm6XtyzTll2N2nuAf+chHMlvHKhfLMU+7DQb7pdOvDUzN5y5HijGoDYVU0tJSAzvSHOvYC2JZYYSa4zDFZO9MlFGWUceaPpc7lavW2vuRvVFyQLt8tp3BWM+2oy8tS5Zq9Gyl+jh3eUB6x3h10+7LGH55PvfS3Ek6dpsKF66b33smkCaQJpAmkCaQJpD2IZA6v3bpkX++S8lUm6+zO7KeDXLBMhva6HDMohZvtIHOcg3Sjk4//H4JoXgejI0nibC2UGqTZlKM+GrplEsZx6v2mg9tEbx0CjLmo22DTtImshsU2sZoGXNud6ax6mlI/aVF2dFYSwfdpNYAkWk6fSe/XMPGEsxbbLQ0umsBxjhaRow0Pq608motjWOLIw1LjVCwiGg5nO1IZI3U8j46I6V1OVIM+iP2Xip5raMJpAmkCaQJpAmkCaR9CKSlS43yha3lkExvIioWe7sdwcbyHfeQQFrObA5ONlldE2ELkNbGkeUw8fnDLGtd9CrfoqyWOg9JBxiW5tkSyeymkx2F3WKp3UzThkHKiH92RMhuCtwMk8aEdcNTAVthRufuLOONvShz4GREL2tJKjtiZjvNsfx+MwJcfixZ+saOjjGBNIE0gTSBNIE0gbQ/gfR6le3qHsHGKDtOoW6a53DLJtq3qHjZZPswsaOwW+y0NmW8HQO7fUg4JKLWWm1nde26dKndLLJ3PWxC/Yu/WXYcrbY7zNrDNc/ZPtN9yFF+N2/cHLB2c4wJpAmkCaQJpAmkCaT/3kDao0UPKfNR9vKfXba48qadDkvSLbPqey9HjP/dJrX3wudeRqu9a3LHrrdwuMfB5ei1N4E0gTSBNIE0gTSB9NsCpMPqYIuiN1/HD5fp7c+3o3Qv3rPbm/0hJ3nXmuzRWV/HsHM7n3t0xyPu98j43PuE/uEG3rtJd/Rj9N69awJpAmkCaQJpAmkCaX8CaS+SH9nc9B7Nvxe1bn8Jfl3stP1Ve4+RwF78+3Cn+/cYzOwd6keM273IvmOe6BEHJIeLliMIeDbVuMcBegJpAmkCaQJpAmkCaf8AaY/vr3vh5v/bTPThvnwfsQccgSPu5sFHP9O6xTWPgP7RQ/cIcgiOYH78yNI7DjdYOkq7TyBNIE0gTSBNIE0g7SsgvV4fel/f2cYdkxH3zvAh0xAPOcH9usdXR4PtvaBx7x8tjozbHZV5uEC6affVZUeJ58MdgPYS2u3RUSeQJpAmkCaQJpAmkPYzkD63UW7aaVXGbts0b+lyWX9ZeW2xx46T47s9WTI81gy3BcwfHyxry6e3q+ymjeUuh+suu1Vb43O5nnkoecsaqj269VoXa6fO7NEFl9TGkqo/OVg2l6PvReq1mmPnpt1W2u/G9pp0S7KbS7wONwTascL2bQgmkCaQJpAmkCaQJpD2D5B2s8Hm/s5ajh2lx/bfa+fm7uVb+1Du0swRH9su73G+Mt6GjduA8lOf+lR7UC4PcvOkfQ+3YGMMDTcdPENliL9UQgpd2mm7bdbItllPOyh1OvVgLK9a7mi5mxOsQXENn0tVjM0oN4ewQ6p3ueVT1D65Ku1x2fC0dgzzbnQG9obFOyixUxU7DPJb3/pWO4BvH0eGgJ0rM3Z1bJvIHGm5M+ZeshwG/D6/yzr5tT3x2599SXwCaQJpAmkCaQJpAmkfAmlAaIQWYw/oDiDpqLxvf/vbf7Eq4xiSDL/bbOw4SDhpo58bjVPfOjYv5pZv/3uxdGfGdJDeDTfccP31139kVT784Q9/6EMfum5VOkxFj5t7iw99NUzg4fOrY3qX+0Fjpgrj2OaxSfdNBw/63dFIY3tElP2JIHcE8o9+9KM3rIqbjsfza47VMS3Lg343bbQM/2Jy+FPnIyJL9g9+8IMEpxm/DkfZLZRdbrjZk7UDF3GbYt2MsWkcj7cGpDF8LI/Wa4vMjiokY4dq/+hHP+r0vh//+Me33347pxpH5e0YrZGlI4LaNR4n7jHj4dhRfRnG7CUuGkNY7rcMsxNwnK3YSZbLcxwnkCaQJpAmkCaQJpD2J5DG7s/jxZ3MnWnR+cck/M53vnPbbbfdeeedP/zhD3/wgx+492QcWztO8lgz/DjDWBkBTKUnOlJzzY8LG7bMoWuSTrMxtHQK4Dj/j5Y7T8WfYUl3BNxx0raRohMKEe84lmEeTzzvbb6TetVcHn684ybUWpGoIwDV9yuC+BmsXnvttfiknI4BRvm73/0u3xpa3TxocLmZeCHQOIS4CMSvFKuLK6644rLLLrv88sv1QoqGqi1hQGXs0z3O+SOajtLwBz7wAVcj1I033kjtw3ZLl13GFa4YU5PyCd7BzOoYL0gKQj/5yU/uueee++6775FHHnnssccefPBBrsVNN4+OHigist5x0tHRDRaedNxjHti50aF37ePHju7UgXxM7EZbpBBE1mCUsCTVVh12GYeUqz/0OYE0gTSBNIE0gTSBtA+B1MG0qLj2+huKOAEfIjnw9Dr7s5/9jAruvvvuO+6445ZbbkERud4pl0HIQBFWmJMNaJNz916ro+CES9XqKAdFJBh/fuMMtjFRjj7KDIzyNddc8573vIfrXH311YRHkIQFM27I0knDrul605M6MLh3X2yr+cFVwao/Pf/+979PfINIJwdjoLf/xp3lxPpgVR0MwB7VoaAC2dn+mlVxg3k6Zx62QR9xccK9995LvQYpwVL+tNTniBBwOE5N1henT3tu0GR+ELrwwgsvuOACaiFFzCxf6werxX55UidDQzKeScpvIojzgIQyanlIk+ANQ4PsGEEQ9JyVjRrDNFp5jjjA8B8qJbixg+AHDhx45plnHnroIbITZC3wHvEwqVmcdO9617uoET/pYRwUlH4asFiQ7zUabmJphIWNdzrNRiCK+MUXX3zRRRe5Go+uvPLKnAGdXKtjsBtMJ5AmkCaQJpAmkCaQ9iGQ8vhSAXTszzEl3YHBjA028ANI4HTXXXelBXoncyAc39GH8GQjCSW++93vJjYuYyUOevnGn+ZsCa6M5+rXJnaX79/LYAZLiGOVkQgJRYgPv0fB+/fPf/5z79+YzIlzvkRbW2XdgcFakdENIgheeuml2MYhfUWNfzM5VWBMtFPmRAPBUuqlM6FJaXTFNT1UnzMVrWGJ5QxDKvAqXXAsKHr00UdFC/qiCswsQ4Wgq0dKo9WCqzEtC7GaIMgncMiNeAMscWLcetixxx2PtyY7Cn5F4Xvf+x6N3b0quGJfP6lJWC7LO5HSKYBRS8NZ083L8S6a+SU28HD++edjg/X9hCA+KSTP6bBq1rz99tufeOIJWOJjI0xajsgNWzgBJAZiemrMbZQA31cQGhYqqwNvqnmymYexzLruwwOeuQdHuuSSS/AMQt1gHqLe+c538jTU6AHP2lKmhrE3gTSBNIE0gTSBNIG0r4BUgISJYhj3TYPqz0tnM9eqsVYm50+UXoqh57oPgX2kXwrviRfry1YFccoVCcAh73GPLIHpmoL82buyLhipD+o56NrLt5pk1hHkkM0bLa+KZz9hEod88fnnn3/qqacAgDf4qfCv5IYxgT6cqTiwfpn8zDPPfMc73oE3luZYXtwffPBBN+iQgq5hrMnrMaW+9nW/b+EkwsADDzxAXpw3W+0n2ksPkAPtfdRXGeQefvjhp59+2tUTplqmI+C8gJOjUCw6fLGhB7ybnBV+kJdpdESx9ONKWMTxoH7fKobsJNKKW4NQk+/ExEZYwqQmAYmwMKxTdUh03333adUc9xJII+BUKJMrn3POOZxScxqAkABGkKIdzqoC96BGCslkOh0z4CM4HBPZwckNBmjv/vvvZ3G6opPx2SPfKLzxMBRtAqmbVOfK5aiL3Sm5tBuk+jcAlldddZWrezcYIE4YmUCaQJpAmkCaQJpA2m9AwpNm2vT6Xk6ha46uMVpaepi9ORzO1KcIb42BgYflVcu5Wq20FW+oqdUjjzzy8ssvP/vss1yHpmDMT+jzcrrzkL96RVaTOnoZHSFNHs9jmqSmTc31rq36vIdLea4Cn0CHSQCA93jYJ/9miluvUnpilJsaVs09Qdj11FNPZQNqFXKwroiF4H7VnTfmc889V/jBVLG3zHdc5jSQiCzYePLJJ40OlI4NNZnfc76ISe5IagwTAQ+a8AlaAiQOjaVl7KF5USUz61cT7PFLqihFlYwgqiFINHYUNFKO7vicrkvdSPZ8qOTRSJWabDDyhMgY8zDfYBEG9RO2f/nLX7766qs0Q8ylgZK9uIWZNCzaxHDZDIijY6jijqV3iD3esSqsqRVUIEu05TR9g3IpEayMsi5cAT5I33bbbZTTfDdt9DGGMwcAPa4NneMjSjzHm0I5LQxrOEYfwIrKeGOJMvkzzrHUVxBdTyBNIE0gTSBNIE0g7R8gxSWlt5LHe63GRPIz2bzCch0M4UZ/OKZ3XNLCBRdccNZZZ1EBryJDk8tFNVEfr5gqY45b/83f/A0DuOEx+EP87LPPRlwrNuOvlMKEBSqe4ydkJjDnoOJUSV91iquSaJm8WVq9cwLO2qRz6ZVk1qpAaKScDn8KBk18w8lFF12k64IZhfn5hDdj3J5++unY1gueuVGLpZunXi6b1wXDcIvnnntOwEa6knH7ri/86CvCiFvCiSv2yOJXqi5lZJRizuLYSGEbrtgVY8zBX3FL8FI//YpDKDpw4IDBS+wBLcsYqWTcv1wV90SmUl1QDqEwRgodFQAzhMoU+/jjj//qV7967bXXuBolCHXYa0SJjcLfWRV/QjiGa6h3xC+88MKTTjqJhnFbooY6zM01sYRbDXFeAsqID8vrpRnDAQWWwJqSUaY01BABJCLTNiJYYjIPWb/Z7TEqpZkWibmnNMMZczTQoEb5P18VSigfhZgF4RjmmcZZOhkpFxNIE0gTSBNIE0gTSPsKSHgiIT1qozbmtNc9JqCI5FdeeaU/myelKeo+77zzqOPkk092c80112g+FpMvp9R1prLYiYKefvrpF154AXOMxC9FVoiDol40pGV1+sl782mnnebX/HUZeKigWm+uKKDsvld5/GMGsLXqRR8zjEQLiFANoLK9jmiqifvlzjtpoBd3N6XtFh6gYLA4b1WogirRR4f/tTyJTlrpPeZ/VcAnIHE7qCioQ43Z1Gc5Gi45oDU2rRgvBRZvTWqXOrAEEnlJ1/S3mm5YFLYZgrpou3CIT6BfugbPe/7558Vpd9555zJAyvULzDBMcKozRpC3XBY1Gf36VQGValK1UUAId++999J5sZlqjQINdjrl8VTNfctIpiUN/cRkZ555JqOjSQlIgYHuCMKLSFdmbYusxsBU3MUti4j8yRyGj4ZRz0mKE/6NAkPHHjoNxxhQh9rHt4SxACltYM/oA0tGn6ygJhMY0/HcerZCzebZ6cpA4PnYMXMCaQJpAmkCaQJpAmn/AKlgBgc0SIAm/nqdFd7waWFDbb6xKup4jhCl99kYhRb/ZPhl0mofsLFOX2QrkmmB0LXXXksXXPPyyy/3Z6kMOPauf+KJJ55wwgnoQ0tvsWNeNZPQEefGG8aCLjYA/rjjjjvllFP81I5+aDIYv6eRAWM6pf1ipLHboDpj38NUQyGIkx0b5SwCmDoFii3EGtv9DDCM9+9yLEoAVbMM3RZmK2q2AFsXZYK2WoYJGMnzdq4J58OZ3KuZXcHJjRjg4osvpkOyu+cxGDNeMLnm6uPWn4888oioCWMxmcbqsY8B9C/Kona28FArnWKMfdmaKhKnDwmELfElD2HHAaRsBCF9M4CoROY2/Bh7rYFHGYeix5deeimPTxY6KTrSYw3TZ3EXZRLEcMAQb33rWwGS1XJuvteyNKRE4AwNb+hQC8hh6a677sJV893LBUh9ofEQbPrnQeS+bdBDs/Ax4x7S+A86ealBoc09J5AmkCaQJpAmkCaQ9huQcmJa5ny69GbP7/1UrkO7qGCiz9UEyxh9wBY13boqOmiGcfnxuDd7fHMmNbHe7pCFH62Z8ZZ/xapwVvydccYZwMA5CIOUXtpaaMRILKq5ht7ptWqbRVJpe+qpp7p6QrMtk+5l94knnvi7v/u7l19+mS6arQ7zA0iKmphsllZfuREU4QqfeiztQMDDlrhqXVDrrttyaACpr+YxzCnHdjNpgz45zUWrAqiNRFjCDymyfWkEhZ0DSExbzIa3tsJs+Xff3VtPzgoEV60ddtAEe2NHTqzaWJE19j8iO9cUZXFNRPpCgAGy0zBg46019shyVmrxa7sLiUzKcl6ubmIgXkRXxG97yja+bN1awcazzz77q1/96qmnnqJtCim51k+aN4qN2LjVawEJPEjHxMcccwzVYalUGA0NAaTAPFsbNDHMdUnkIccDgFJsB9obBKmRtimwjYdId+Gq0B6JWBxBTltGhTHR6OCh+sYv17GF6ATSBNIE0gTSBNIE0r4CUos9dEDjDz/88JNPPsldYKb38o4h6aM7FXjbvueee7ytPv74467u204oIC33+mlfmzZe7Es5hpqYZl2ckR83fScmD22SpJ0TWRQRPbbqaUyt9mrLhJoAW9m0hoA++bfNDXFUKw1Uc+x5F3/11VcpvanetXNWFE+SnQ0oS8PCFZ7N8xjbqzxh2QPbnpRqiedmopVWT429aVilEMhNX/rLqcUqAXHOxoyhjvotltacN5x11llCRwYe2z4ut7xp9cs4vGRkdrZTgL5wEhrPOeccyumzvQqFu20WWYyU62jlScm46vTxw9X96aeffu6556rDUw8cOMBkJQrTcyMj3fbxYOwYXnzYVkQlmtTjWO4FBlyIIZ555hnaxpvxiEJYrSGvlBodJeZYOMS1DNZMgwLeaN7zBx54AGPcidsYQzmSru9fFX1xJ/xTCP49QTyN9X2i/FQjmobU1RaT7BUzzQgQWaftptScvtjPc8NHI+A4q2YCaQJpAmkCaQJpAmn/AKk5VjfMDCRQVIJju9j4VQda6ok8APPoo4++8sorr7322vPPP++d1RM1W3neO+iIZwjTLsn4KPwQovSx38uxGyooZ7GlRx1tUmnXG2SX2/ClgmYkgxwzl9nQihru1YR7J/a1Gke/eGY81AoJlukXGUx3dPTXf/3Xf/u3f+slG6k+CfRRXPMXX3wRkFCgPuZ5+9vffvbZZ/OqqDXWFM+MQKiQj2jEbNtHZFXjByDE6uqUzUnt6tBA664ZqZ2VstNyaXQ7TmrSUnA35bxSo47aUwkORV966auDTgt48FOOwtggOxR5ksdjgLbZyw0ioM7PMGZwMRLhvzinYYUTt0NTx9ctnd4NmsaIS1cFY208xBdzA/KyV8OKAcvVr63/CUjLIamRoggWDzptm0sUWmEluDr++OPf/OY3k5EqylFQgRppA4Z5WnPxY0l8wxANG7kE5DSjx5brF103gvdvoENiNA+Q2AgRhs5xVuIE0gTSBNIE0gTSBNL+AVKlXEPCeHN96KGHSkGgDlpoorl0RooTPnkxVefOO+9seXbRy9jRbwmkzkFh7PbTYZKXX36Zv4bVZiSb407FCsdtFdA4oHcAaZzt0dkwOvIWixPgf+yxxziWTknYvC2PjGdEWAjYtG1t0tgrcICzUOqFF1745S9/CTbFb8T3EzE9AX7mR4FvgYcIDaJaH1VkOM6uK3HB2zYe4A1gaImwbbuD+Y4SxEbv/W3W495zDBd3tUf5OLp4hJ26M3KhxgQdgoI4MPRdH1kCYps+abLNj7iFqKnV0a2NT/a+QKhc2ND25W2zWMikIMh8XMLgwmtbf04Q4yCd8Ba/8rzlKSyIlKFisDh5VUjEWE06twNU6QvNwrfXEuLAWZLoOMRx0EwPZOfKhHJtH3ZP6OqMM8445phjTjrpJONa+QfoEIdOeDy7o6yyjpbr9t03303n2BO54cFg2uomKqU0tqaK/oUQlvfyNDU1bKQbmxxNIE0gTSBNIE0gTSDtHyAtd37hl81l046qOQ0WvYOeeuqpXihbcNseQzrrdN6+Q6/tfJ13lgyBs3Jh2zqHIrDrSV/9XdvIpm/wLdRpqjQ3HSc0d6Zge0CXakBHz6wK6yLCUZgNqxh+29vexpaedGgKOm0zPYC0XHXs1zYM4qPY8JD4NFi2YnsmtgS6o/KI70/maXHUWrauVlwQkI499tizzjqL0iCzfRtz3z42tMInX0STzcoPbvI3kcf+mAGpGKlc2HaQZAVNdOcGfUrmas8995zRyhXSSM1wpRev7czen+0e7lrOCmqNRBdeeCE/C/Aq9OmCILTEjoxY+i/tjdijWXUhmQgNDo0jrMCF2nEphHTiIBW17ArDBZAcoBB3bOI5AqSSMLgoz3xwVYzmem/LAxbXHWNpQjn03N6j5PWk83Ioebm5wDCTVqzjV0ozVgIhYccgeNpppwmfeAKyfv3FL34BSPyEBUee7v871mUCaQJpAmkCaQJpAmkfAGlttrrUUjd9gOeI/OCUU07xpnv66afTCEt79WS2Xso7TKV9ZMac8vKA2yYQizfq2LWNC8eakxaHtLlPM5vhsHnwEXotj4smhhv1ORO79rrcnoBXXHHFRRddxAla5oRyH6cRKbIah/aN6dr8oANj2lfR2zY6xhFxi/qdv0IzRW7tUc78Y2/H5VkpORMLcccTTzyRJxUk5A1jvdO4p4QWkJeW0aTq5rl6XzhYYrgp5hY1EdOf2DOaBKFXXnlFVMPtME8PI6dhRHEj4Gzpu7Ytr0KzdWLFSNRIina7FjMYhjQfG0qqP9JERrTQUiuyaFUCKJ/p80a5GqWSlLbSOYiNqprTQ3kqY9J/nDkU5rFKqDaIb0ekTivsYBujkgpYQqHtSv3Uqv6RTr15+rg6bRZ///33q9nCNt7e/lbsyPpwLn5+8cUX3XT6C8UulTmBNIE0gTSBNIE0gbRPgLTce7Fvum2sjBUI8erZPizvWhU3uZfnIOT19LrrrqPc9oVcO414HPjXpDPNUiiLMsxll13Guh2e0QLslNtSpdSdMy2NNF7r8dZBzmoGV9VaUtUJNGyGTsmgAh52clMiaWuDl8d7fHFVWuqjLb8UTnQQL3nbGTBSMMOKLUoRShXpZaTl4SvxiQcq8qrdAOQ+j2/dl2sx2A2r0k49fi06WpN6bY/IsgG4u5EuHeqO1do8UfSCt75MIMhqFM5kLXBaHlVYocAOg+nDAx9oMp2hg5OR9Mwzz+RV6OiuRWuaFFMtT7JZ7kjV2TPwQ7p2mBqqTrriPVgq3mY+HJY5OgblNbsXYrUmzcjYbt3s21R4q90wUPqy3im5eJ4/D1bXsFQqdmvwXFHosBn/Kty0SB5oRVAGptIaOpNmTZkTSBNIE0gTSBNIE0j7DUgt6O10k6Ka/DWXzeqY+9DB4r4j3CLat+cd3+nJT7Cx3w3bMEmxlucFS63n6drs8HLr57XjPVTQL37yy7BXkxZBtQOj9+87ViUV54XLLSyX63xaA5ODBvjO8SVjx4TcfPPNRUfUCkveqjuIur2TlnO149TkjgBpnVXH+yVskYBfP7wqQYj7tmB7eQDeGpDKmSigask3z9ZWd2MQxCcKgaFdFFVjrAKPpdMPILXZJSBxx7DUvLBBpC11EKEQFmwr9sLXEpRv2iilBXcydxK1p1LR8nIN1ThPaBwvTbS1kW7N7gUn7WvZ0dEswhCdot2xNy3Fp5YG0zVlrp2Z3ZnKKEBR+7Cjg8/ixvJrBeE//elPGb1EnI4JX4P6BNIE0gTSBNIE0gTSPgHS8v27qCZC46W/l9pcrfikT/uJ1JqZAaE1wydAO4CzK6t48878nX9GgwNFERwzpGuwXH6PR1ATvtjxz+1JVHJGi5xLI0jFY6fv8imXrC5lbyulXu6b3HfNv1uInqcyWIegtFlPuxEtUyGX/jSMWjpF3+wbLwjuYWNBn/zHAvJWvyypLYl3DoqGItixeVAu24eBTtSmYXo2Ioij/Fkyx6Y+R1ILBBZvtP06dTWyoFbogtu8vBBuk8m1aeWx/1FGKRbq28ZY4Y/hPoREf8yPj5FubbPR8NmXD5Qp/6ur0ljW0USN+67jhKEdy+b5gjDZoY8s267fd64KLPW8uYMsu/y3EbUJpAmkCaQJpAmkCaR9CKQRKfW5On0RXnvXcd7wOAolcmsU1xy0jknFHk3yhp/eoZuPLsIZp8KscbmjncpbbZl0CY4NAdTakua+c3ds826zn0s+B1l2KiZsXj6jppbWIXvXb01yu08uUbQWJo2tJ4cTNwqMvXsG/8ttJddork1YRyH/o0lX/pd7NbgsQ1naHkHscsZ/c/Rsl59bVoXSOltl7LU0FlSPuGJT8DUbVfqskt6ywkhTSA/ji0Wq2G1QXtp9hNCbPzVUjeF4N0da89Im4knddpDdtCKuM3vGgvOlNtY4nECaQJpAmkCaQJpA2p9AGnPWy51rlgbeMR7YwuX4LJ33tBAlIZdljbntZSw4T++DSTZrZVFnzOTra1mVuwFpmcUYvJdTnEs97Ojuy6hjXOMzfxqhZg+XE/Frn/x3hOWm7OlzDEN5WIiqLFW9FnBujp7NUN+8KqGo8XQp7+cPVTaxNJZRRartyxM/LMX2Fr/f/PixDKHHmN7zMRYvkb9p7k1VjAP8miAYpaVWib800CaTE0gTSBNIE0gTSBNI+xNImxpZ23h6u5k3TbU0/1LIPWJmM0jYrPCF3yxfXJS1cWHvfC6nODdHmR3dfbeoJmdaDkObut0cpLYreanPpd/scWzaHD1HPJzrLFW3G9t74XMtzWWMLztyu1vbtRSHZbS29nBpr+0etcne2n+L3cqOg+YE0gTSBNIE0gTSBNJvEZC26G57yLGj+fdi1x29c1OALQ23v8dvIbj2fn/Evez2Ln6UfG5X5iEdaHsZY9BewuDDwtJmfHLI0fPInG0vLrSb7FsGsr0YdwJpAmkCaQJpAmkCaX8A6f8CGyNfgXju8cYAAAAASUVORK5CYII=",
      "text/plain": [
       "BufferedImage@105a53e6: type = 10 ColorModel: #pixelBits = 8 numComponents = 1 color space = java.awt.color.ICC_ColorSpace@7d546d71 transparency = 1 has alpha = false isAlphaPre = false ByteInterleavedRaster: width = 840 height = 168 #numDataElements 1 dataOff[0] = 0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val onerow = svecs.view(28,28*opts.dim);\n",
    "val nc = onerow.ncols;\n",
    "val tworows = onerow(?,0->(nc/2)) on onerow(?,(nc/2)->nc)\n",
    "show((tworows.t*500+128) ⊗ ones(3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting Right Singular Vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we obtained the singular values and left singular vectors from the model's modelmats array. The right singular vectors grow in size with the dataset, and in general wont fit in memory. But we can still compute them by running a predictor on the dataset. This predictor takes a parametrized input file name for the matrix $M$ and a parametrized output file name to hold the right singular vectors. A key option to set is <code>ofcols</code> the number of samples per output file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    },
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val (pp, popts) = SVD.predictor(nn.model, dir+\"data%02d.fmat.lz4\", dir+\"rSingVectors%02d.fmat.lz4\")\n",
    "popts.ofcols = 100000                  // number of columns per file, here the same as the input files\n",
    "popts.nend = 10                        // Number of input files to process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicting\n",
      " 2.00%, ll=-0.01308, gf=1.979, secs=0.6, GB=0.06, MB/s=98.93, GPUmem=0.99\n",
      " 3.00%, ll=-0.01367, gf=2.877, secs=0.7, GB=0.09, MB/s=143.85, GPUmem=0.99\n",
      " 5.00%, ll=-0.01273, gf=4.519, secs=0.7, GB=0.16, MB/s=225.94, GPUmem=0.99\n",
      " 6.00%, ll=-0.01348, gf=5.316, secs=0.7, GB=0.19, MB/s=265.76, GPUmem=0.99\n",
      " 8.00%, ll=-0.00652, gf=6.551, secs=0.8, GB=0.25, MB/s=327.52, GPUmem=0.99\n",
      "10.00%, ll=-0.00664, gf=7.503, secs=0.8, GB=0.31, MB/s=375.12, GPUmem=0.99\n",
      "11.00%, ll=-0.00559, gf=8.061, secs=0.9, GB=0.34, MB/s=402.99, GPUmem=0.99\n",
      "12.00%, ll=-0.00574, gf=8.593, secs=0.9, GB=0.38, MB/s=429.59, GPUmem=0.99\n",
      "14.00%, ll=-0.00634, gf=9.483, secs=0.9, GB=0.44, MB/s=474.13, GPUmem=0.99\n",
      "16.00%, ll=-0.00670, gf=10.283, secs=1.0, GB=0.50, MB/s=514.10, GPUmem=0.99\n",
      "17.00%, ll=-0.00566, gf=11.113, secs=1.0, GB=0.56, MB/s=555.59, GPUmem=0.99\n",
      "19.00%, ll=-0.00624, gf=11.394, secs=1.0, GB=0.60, MB/s=569.64, GPUmem=0.99\n",
      "20.00%, ll=-0.00687, gf=9.531, secs=1.4, GB=0.66, MB/s=476.53, GPUmem=0.99\n",
      "22.00%, ll=-0.00606, gf=9.794, secs=1.4, GB=0.69, MB/s=489.65, GPUmem=0.99\n",
      "24.00%, ll=-0.00577, gf=10.186, secs=1.5, GB=0.75, MB/s=509.23, GPUmem=0.99\n",
      "25.00%, ll=-0.00639, gf=10.744, secs=1.5, GB=0.82, MB/s=537.13, GPUmem=0.99\n",
      "27.00%, ll=-0.00755, gf=11.026, secs=1.5, GB=0.85, MB/s=551.25, GPUmem=0.99\n",
      "28.00%, ll=-0.00601, gf=11.259, secs=1.6, GB=0.88, MB/s=562.87, GPUmem=0.99\n",
      "29.00%, ll=-0.00534, gf=11.369, secs=1.6, GB=0.91, MB/s=568.40, GPUmem=0.99\n",
      "30.00%, ll=-0.00536, gf=11.405, secs=1.6, GB=0.94, MB/s=570.18, GPUmem=0.99\n",
      "31.00%, ll=-0.00609, gf=11.506, secs=1.7, GB=0.97, MB/s=575.24, GPUmem=0.99\n",
      "32.00%, ll=-0.00692, gf=11.958, secs=1.7, GB=1.03, MB/s=597.85, GPUmem=0.99\n",
      "34.00%, ll=-0.00610, gf=12.029, secs=1.8, GB=1.07, MB/s=601.38, GPUmem=0.99\n",
      "35.00%, ll=-0.00552, gf=12.483, secs=1.8, GB=1.13, MB/s=624.08, GPUmem=0.99\n",
      "37.00%, ll=-0.00629, gf=12.579, secs=1.8, GB=1.16, MB/s=628.90, GPUmem=0.99\n",
      "39.00%, ll=-0.00694, gf=12.909, secs=1.9, GB=1.22, MB/s=645.40, GPUmem=0.99\n",
      "40.00%, ll=-0.00606, gf=12.900, secs=1.9, GB=1.25, MB/s=644.94, GPUmem=0.99\n",
      "41.00%, ll=-0.00533, gf=11.435, secs=2.2, GB=1.29, MB/s=571.70, GPUmem=0.99\n",
      "43.00%, ll=-0.00585, gf=11.732, secs=2.3, GB=1.35, MB/s=586.55, GPUmem=0.99\n",
      "44.00%, ll=-0.00650, gf=11.901, secs=2.3, GB=1.38, MB/s=595.02, GPUmem=0.99\n",
      "45.00%, ll=-0.00678, gf=12.242, secs=2.4, GB=1.44, MB/s=612.03, GPUmem=0.99\n",
      "46.00%, ll=-0.00566, gf=12.330, secs=2.4, GB=1.47, MB/s=616.44, GPUmem=0.99\n",
      "48.00%, ll=-0.00544, gf=12.436, secs=2.4, GB=1.51, MB/s=621.76, GPUmem=0.99\n",
      "49.00%, ll=-0.00614, gf=12.633, secs=2.4, GB=1.54, MB/s=631.58, GPUmem=0.99\n",
      "50.00%, ll=-0.00663, gf=12.626, secs=2.5, GB=1.57, MB/s=631.24, GPUmem=0.99\n",
      "51.00%, ll=-0.00671, gf=12.716, secs=2.6, GB=1.63, MB/s=635.76, GPUmem=0.99\n",
      "53.00%, ll=-0.00531, gf=12.772, secs=2.6, GB=1.66, MB/s=638.52, GPUmem=0.99\n",
      "54.00%, ll=-0.00554, gf=12.914, secs=2.6, GB=1.69, MB/s=645.61, GPUmem=0.99\n",
      "55.00%, ll=-0.00623, gf=13.004, secs=2.7, GB=1.72, MB/s=650.13, GPUmem=0.99\n",
      "57.00%, ll=-0.00698, gf=13.277, secs=2.7, GB=1.79, MB/s=663.77, GPUmem=0.99\n",
      "58.00%, ll=-0.00614, gf=13.361, secs=2.7, GB=1.82, MB/s=667.97, GPUmem=0.99\n",
      "59.00%, ll=-0.00569, gf=13.541, secs=2.7, GB=1.85, MB/s=677.00, GPUmem=0.99\n",
      "60.00%, ll=-0.00531, gf=13.465, secs=2.8, GB=1.88, MB/s=673.20, GPUmem=0.99\n",
      "61.00%, ll=-0.00623, gf=12.464, secs=3.1, GB=1.91, MB/s=623.11, GPUmem=0.99\n",
      "63.00%, ll=-0.00680, gf=12.573, secs=3.1, GB=1.98, MB/s=628.60, GPUmem=0.99\n",
      "64.00%, ll=-0.00570, gf=12.741, secs=3.2, GB=2.04, MB/s=637.00, GPUmem=0.99\n",
      "65.00%, ll=-0.00590, gf=12.774, secs=3.2, GB=2.07, MB/s=638.62, GPUmem=0.99\n",
      "66.00%, ll=-0.00621, gf=12.888, secs=3.3, GB=2.10, MB/s=644.32, GPUmem=0.99\n",
      "68.00%, ll=-0.00668, gf=12.961, secs=3.3, GB=2.13, MB/s=647.97, GPUmem=0.99\n",
      "70.00%, ll=-0.00667, gf=13.127, secs=3.3, GB=2.20, MB/s=656.26, GPUmem=0.99\n",
      "71.00%, ll=-0.00588, gf=13.283, secs=3.4, GB=2.26, MB/s=664.09, GPUmem=0.99\n",
      "73.00%, ll=-0.00620, gf=13.354, secs=3.4, GB=2.29, MB/s=667.62, GPUmem=0.99\n",
      "74.00%, ll=-0.00677, gf=13.412, secs=3.5, GB=2.32, MB/s=670.51, GPUmem=0.99\n",
      "76.00%, ll=-0.00665, gf=13.559, secs=3.5, GB=2.38, MB/s=677.86, GPUmem=0.99\n",
      "78.00%, ll=-0.00564, gf=13.674, secs=3.6, GB=2.45, MB/s=683.64, GPUmem=0.99\n",
      "79.00%, ll=-0.00625, gf=13.753, secs=3.6, GB=2.48, MB/s=687.60, GPUmem=0.99\n",
      "80.00%, ll=-0.00663, gf=13.699, secs=3.7, GB=2.51, MB/s=684.90, GPUmem=0.99\n",
      "81.00%, ll=-0.00736, gf=12.811, secs=4.0, GB=2.54, MB/s=640.48, GPUmem=0.99\n",
      "83.00%, ll=-0.00581, gf=12.928, secs=4.0, GB=2.60, MB/s=646.36, GPUmem=0.99\n",
      "85.00%, ll=-0.00618, gf=13.078, secs=4.1, GB=2.67, MB/s=653.81, GPUmem=0.99\n",
      "86.00%, ll=-0.00659, gf=13.167, secs=4.1, GB=2.70, MB/s=658.28, GPUmem=0.99\n",
      "87.00%, ll=-0.00715, gf=13.275, secs=4.1, GB=2.73, MB/s=663.66, GPUmem=0.99\n",
      "89.00%, ll=-0.00580, gf=13.420, secs=4.2, GB=2.79, MB/s=670.92, GPUmem=0.99\n",
      "91.00%, ll=-0.00605, gf=13.475, secs=4.2, GB=2.85, MB/s=673.69, GPUmem=0.99\n",
      "93.00%, ll=-0.00704, gf=13.566, secs=4.3, GB=2.92, MB/s=678.25, GPUmem=0.99\n",
      "95.00%, ll=-0.00573, gf=13.655, secs=4.4, GB=2.98, MB/s=682.68, GPUmem=0.99\n",
      "96.00%, ll=-0.00579, gf=13.692, secs=4.4, GB=3.01, MB/s=684.53, GPUmem=0.99\n",
      "98.00%, ll=-0.00652, gf=13.777, secs=4.5, GB=3.07, MB/s=688.77, GPUmem=0.99\n",
      "100.00%, ll=-0.00664, gf=13.789, secs=4.5, GB=3.14, MB/s=689.38, GPUmem=0.99\n",
      "Time=4.5490 secs, gflops=13.79\n"
     ]
    }
   ],
   "source": [
    "pp.predict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "BIDMach's actual SVD model is only slightly different from the code presented here: First it computes the subspace updates on *minibatches* of data for the first few iterations to more rapidly get near a good model. Secondly it alternates subspace iteration with Rayleigh-Ritz iterations for faster convergence."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "name": "scala",
   "version": "2.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
